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FOREWORD 

This is the first of three volumes describing the Performance 
Analysis and Design Synthesis (PADS) computer program. This 
volume is devoted to a complete program formulation. Volume II 
contains programming and numerical techniques and Volume III 
is a user manual. 

The development of PADS was conducted by McDonnell Douglas 
A stronautic s Company at Huntington Beach, California, under NASA 

Contract NAS9- 12059, under the cognizance of Mr. Robert Abel, 
NASA, MSC, Houston, Texas. The key MDAC personnel who 
formulated and programmed PADS are Messrs. Murray H. Rosenbe 
John W. Hensley, and Michael Beach. Valuable programming 
assistance was given by Larry Ong, Fred Gangloff, and 
Sheldon Herman. 
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ABSTRACT 

The Performance Analysis and Design Synthesis (PADS) computer 
program has a two-fold purpose. It can size launch vehicles in 
conjunction with calculus -of-variations optimal trajectories and 
can also be used as a general-purpose branched trajectory opti- 
mization program. In the former use, it has the Space Shuttle 
Synthesis Program as well as a simplified stage weight module 
for optimally sizing manned recoverable launch vehicles. For 
trajectory optimization alone or with sizing, PADS has two tra- 
jectory modules. The first trajectory module uses the method of 
steepest descent; the second employs the method of quasi- 
linearization, which requires a starting solution from the first 
trajectory module. 
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Symbols 

A 

A 

e 

A 

exit 

a 

v 

a 

V 

a 

a 

m 

a 


a 



Description 

A matrix (steepest decent formulation) 

Net nozzle area 
Nozzle area per engine 
Acceleration vector 
Component of a in v equation 
Component of a in Y equation 
Component of a in ^ equation 
Component of a in m equation 
Semi-major axis 

Coefficient in booster- stage empty-weight 
equation 

B matrix 

Coefficients in orbiter-stage empty-weight 
equation 

Drag coefficient 

Zero-lift drag coefficient 

Lift coefficient 

Uncorrected lift coefficient 

Lift coefficient at a = 0 

Lift coefficient slope 

Moment coefficient at a - 0 

Moment coefficient slope 

Constant lift term in governing equation 
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Unit 



ft 


/ rad 

/rad 

lb 
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Symbols 

Description 

Unit 

C 

a 

Constant angle of attack 

rad 

c. 

1 

th 

Multiplier for 1 homogeneous solution to 
linearized state and co- state equations 


D 

Aerodynamic drag 

lb 

D ( ) 

Direction cosine matrix 


D' 

Diagonal of A matrix 


D b 

Base drag 

lb 

D .. 

y 

Partial derivative operator including 
effect of steering vector u 


D 5 * ' • 

Partial derivative operator, steering vector 


y 

u fixed 


d Re£ 

Aerodynamic reference length 

ft 

(dP ) 2 

Metric of control and parameter changes 


E 

Energy or eccentric anomaly 

ft 2 /sec 
or rad 

e H 

Unknown control functional 


e r 

Earth radius 

ft 

e 

Orbit eccentricity 


F 

Vector comprised of state and co-state 
differential equations in QL module 


f rated 

Rated vacuum thrust (total) 

lb 

F 

VAC 

Vacuum thrust per engine 

lb 

f 

Differential equations of motion 


f 

Nominal differential equations of motion 


G 

Diagonal matrix of coefficients that 
multiplies acceleration vector 

xii 
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Symbols 

GM 

G. 

J 

g 

g 

°max 

g r 

H 


h. 

i 


sp 


sp 


B 


3 Peff 


S P, 


V 

V 

J 

j 

K 

I 

K 

k 


Description 

Gravitational constant times Earth mass 

Variable part of cut-off function for arc j 

Gravitational acceleration 

Maximum total acceleration 

Reference gravitational acceleration 

Momentum (ft /sec), or Hamiltonian, or 
matrix of homogeneous solutions to linear- 
ized state and co- state equations 

Altitude 

til 

i c n homogeneous solution to linearized 
state and co- state equations 

Vacuum specific impulse 

Booster vacuum I 

sp 

Effective I (with throttling) 
sp b 

Orbiter vacuum I 

sp 

Inclination angle 
Unit vector along latitude 
Unit vector along longitude 
Functional to be minimized 
Control blend factor 

Vector of governing equations with elements 
K r k 2- k 3 

Ratio of dual engine I 's 

Induced drag coefficient 

Expression used in vertical rise or 
pitchover control governing equation 


Unit 


■ 2 
ft/ sec 

ft/ sec^ 


ft 
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Symbols 


Description 


Unit 




L 

L 


u 


max 


M 


*^CG 

m 

N 


N, 


Nn 


N. 


PL 

P 


Pa 

Q 

• 

Q 

q mult 

q 

R 


Expression used in vertical rise or 
pitchover control governing equation 


Aerodynamic lift lb 

Uncorrected lift lb 

Maximum lift lb 


Mach number, or minor of A matrix, 

or Lagrange multipliers for state boundary 

conditions 


Aerodynamic moment 

Aerodynamic moment about center of gravity 
Mass 

Minor vector of A matrix 

Last subarc of trunk 

Last subarc of first branch 

Last subarc of second branch (N^ = if 
no branches) 

Payload 

Adjustable parameter, or propulsion vector, 
or particular solution to linearized state and 
co- state equations 

Atmospheric pressure 

Heating 

Heating rate 

Heating flag 

2 

Dynamic pressure (lb/ft ) or control 
variables that are not free to be optimized 

Radial distance to vehicle 


ft -lb 
ft-lb 
slugs 


lb 


lb/ft 2 

Btu/ft 2 

Btu/ft 2 /sec 


ft 
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Symbols 


R 

R 

c 

R 

e 

R 


ey 


D 

SFC 

S Ref 


2^ 


T 

t mult 

Tol 

T VAC 


t 

t 

e 



U 


u 

V 


Description 

Ratio of throttled dual engine thrust 
Apogee radius 
Reynolds number (unit) 

Perigee radius 

Cross-range 

Downrange 

Specific fuel consumption 
Aerodynamic reference area 
Total range 

Matrix of adjustable parameter sensitivities 


Unit 


ft 

ft' 

ft 

ft 


ft 

lb/lb/hr 

ft 2 

ft 


Total solution to linearized state and co- state 
equations 

Thrust 

Number of engines 

Vector of constraint tolerances 

Total vacuum thrust 

Time (usually arc time) 

Elapsed time 
Trajectory end time 
Initial time 
Denotes transpose 
Control vector 
Steering vector 
Relative velocity 


lb 

lb 

sec 

sec 

sec 

sec 


ft/ sec 
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Symbols 

Description 

Unit 

V A 

Velocity at apogee 

ft/ sec 

V I 

Inertial velocity 

ft/sec 

V 

p 

Velocity at perigee 

ft/ sec 

r 

w 

Vehicle weight (lb), or control weighting 
matrix, or approximate solution to w 


W BO 

Booster burn-out weight 

lb 

W 

Booster-stage empty weight 

lb 

e B 

W 

e 

Orbiter- stage empty weight 

lb 

w l!o. 

Lift-off weight 

lb 

w 

o 

Orbiter gross weight 

lb 

W 

Booster propellant weight 

lb 


lb 

w 

Orbiter propellant weight 

Po 


w 

In-plane control vector (subset of U) 


X 

Asymptotic injection parameter 

ft 

X CG 

Body x station instantaneous center-of- 
gravity 

ft 

X CGR 

Body x station reference center-of-gravity 

ft 

X E 

Body x station of engine thrust centroid 

ft 

X T 

Tail length body station 

ft 

X 

Body station coordinate along axis (ft) 
or quasitime 


Y 

Asymptotic injection parameter, or vector 
comprised of the QL module state and 



co-state vectors 


M 

Diagonal matrix of parameter weighting factors 



State vector 
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Section 1 
INTRODUCTION 


The Performance Analysis and Design Synthesis (PADS) computer program 
provides the capability to synthesize launch vehicle design in conjunction with 
an optimally shaped and staged trajectory. It also permits generalized trajec- 
tory optimization including branching for a wide variety of endo-atmospheric 
aerospace vehicles. 

The synthesis capability, derived from the Space Shuttle Synthesis Program 
(SSSP), Reference 1, is oriented toward manned reusable launch vehicles. 
However, a simplified synthesis module is available in the program which 
provides a general two-stage launch vehicle design capability. 

The trajectory and staging optimization portion of the program employs a 
closed-loop steepest descent method for an approximate solution. Moreover, 
the steepest descent solution may then be used as a starting guess for the 
program's quasi-linearization algorithm which determines the exact solution 
of the calculus of variations multipoint boundary value problem. 

This document is the first of three volumes, and is devoted to the formulation 
of PADS. Volume II is the programmers document and Volume III is the 
User Manual. The bulk of Volume I is devoted to the development and dis- 
cussion of the trajectory formulation. The first five sections describe the 
types of simulations that are available in the program, deferring the discus- 
sion of control and parameter optimization until later sections. Section 6 
describes the general two-stage launch vehicle synthesis model (also called 
the Phase I sizing module). Section 7 is a brief discourse on the interplay 
between. the trajectory model and the Space Shuttle Synthesis Program. The 
detailed documentation on this synthesis model is available in Reference 1. 

Section 8 describes the auxiliary print computations that are available in the 
program. 
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In Sections 9 and 10 the general solution for the control vector is presented for 
non-optimal situations. This leads to the key variational aspects of the pro- 
gram, where optimal control and staging are formulated. Sections 11 through 
15 show how the steepest-descent formulation calculates optimal steering and 
parameters. These sections include discussions of the adjoint differential 
equations, influence functions and numerical analysis and solution converg- 
ence techniques. Section 16 presents the details of the necessary conditions 
for the exact solution of the multi-point boundary value problem in prepara- 
tion for a presentation in Section 17 of how the method of quasi-linearization 
satisfies these conditions. 
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Section 2 

TRAJECTORY SIMULATION 

In this section of the PADS formulation document, the earth and vehicle 
coordinate systems are defined and the equations of motion are given. In 
addition, the concept of the control vector, U, is defined. 

2. 1 EARTH-RELATIVE COORDINATE SYSTEM 

The Earth -relative flight-path coordinate system is illustrated in Figure 2-1 
below. The earth model is spherical- rotating with a central-force gravita- 
tional field. The origin of the R vector is geocentric. 

2. 2 VEHICLE COORDINATE SYSTEM 

The vehicle is treated as a point-mass moving in three degrees of freedom. 
The applied load directions and control angles relative to the basic coordin- 
ate system are shown in Figure 2-2. It should be noted that the vehicle 
banks around the velocity vector, hence there is no yaw angle. 

“ — CR155-I 



Figure 2-1. Earth- Relative Flight-Path Coordinate System 
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2. 3 EQUATIONS OF MOTION 

The equations of motion in relative coordinates are derived in many other 
sources. They will be presented here in engineering notation. The meaning 
of various terms is given in Figures 2-1 and 2-2 and the acceleration compo- 
nents depend on the simulation model. 


V = Ruj cos p (cos p sin Y - sin p cos ^ cos Y ) - g sin Y 


+ — 
m 


(2. 3-1) 


Y = 


Rgj 

= w cos p 2 s in 4 m — y- (cos p cos y + sin p cos ^ siny) 


+ cos Y 


( v - g ) _£L 

R. + V 


(2. 3-2) 




cocos 

cosy 


Rw sin P sin 4 1 
V 


- 2 cos ^ sinY 


+ s in p 


V cos Y sin 4* _ \ 

^ + 2ojJ + 

R cos p / 




V cosY 


(2. 3-3) 


R - h = V sin y 


(2. 3-4) 
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P = COS y COS 

XV 


(2. 3-5) 


P- 


V cos p s in ^ 


V cos p 

The rate of change of mass is 


m = a 


m 


(2. 3-6) 


(2. 3-7) 


The equation for heating rate is 


Q = Q 


MULT 


17600. V7 1 (j^) 


3. 15 


(2. 3-8) 


Particular terms in the above equations are defined below by equations or 
functional dependencies. 


GM 

g = id 


(2. 3-9) 
(EQUA3, STATEF) 


The following functional dependencies are characteristic of the type of aero- 
dynamic and propulsion simulations. * 


(2. 3-10) 
(VT, UT) 

(2. 3-11) 
(VT, UT) 

(2. 3-12) 
(EQUA3, STATEF) 


LIFT: L = 

<! S Ref C L 

ja, M 

DRAG: D = 

<l S Ref C D 

I a, M j 


BASE DRAG D, = D, { h) 


We have purposely omitted the functional dependencies of T, 6^,, a and 0 and a, 
which will be discussed in Section 2. 4. 


*T, L, and D are model-dependent force terms described in Section 3. 


2-3 



2. 4 INTRODUCTION TO THE CONTROL VECTOR 

If we examine the equations of motion we find that they contain kinematic and 
dynamic terms. The dynamic terms involve the applied loads. In a general 
sense, the equations of motion may be expressed as 


y = f { y I + G [ a ! 


(2.4-1) 


where, a, called the acceleration vector, may be defined for example as : 


a = 



T cos (a - OjO-D - cos a 


m 


T sin (a - 6g) + L - sin a 


m 


T sin (a - 5 E ) + L - sin a 


m 


m 


cos $ 


sin 


(2. 4-2) 
(ACCEL, 
APPLY 
VT, UT) 


and G is a diagonal matrix whose elements are 


[°] 


t 0 


0 


1 


v cos \ 

o o 


o 

0 

0 

1 


(2.4-3) 


The functional dependency of "a" may be expressed 


a = a { T ’ 6 E’ a ’ ^ ’ y ’ t | 


(2.4-4) 


'■'The example given is for the single -engine moment balance simulation. 
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Of the dependencies exhibited in Equation (2. 4-4), it is apparent that the 
"decision quantities" are T, 6^,, a, and 

We group these into the vector U, called the control vector 

U T = (T, 6 e , a, 4>) (2.4-5) 


and using this notation rewrite Equation (2. 4-1) 


y = f {y [ + [G] a ju, y, t} 


(2.4-6) 
(DER3A, NEDRV) 


For convenience, in later discussion, the control vector may be divided into 
subsets. We class T, and ci as the in-plane control vector, w, and define 

hj 

the steering vector as: 


u T = (a, 4> ) 

As has been described in Section 1, the purpose of the trajectory optimization 
portions of PADS is to solve for the time history of U that satisfies all alge- 
braic and variational problem constraints. 
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Section 3 

APPLIED LOADS MODELS 

This section describes the various applied loads models currently available 
in PADS. The applied loads are divided into aerodynamic and propulsion 
models. 

3. 1 AERODYNAMIC MODELS 

There are three types of aerodynamic models available in PADS, These are: 
(1) asymmetric linear lift variation with quadratic drag polar; (2) bivariate 
tabular lift and drag coefficients as functions of angle of attack and Mach 
number; and (3) static moment balance aerodynamics, 

3 l # l Asymmetric Linear Aerodynamic Model 

The equations for total lift and drag coefficients depend on coefficients that 
are input tabular functions of Mach number. 

(3. 1-1) 

(BER0CO, 
AER0CO) 

(3. 1-2) 
(BER0C0, 
AER0C0) 


D 


D 


+ k C. 


O 


= C, 


a + C, 


J o 


The aerodynamic lift and drag are calculated from these coefficients as: 


L - C l q S REF 


(3. 1-3) 
(VT, UT) 


D - C D q S REF 


(3. 1-4) 
(VT, UT) 
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where 

(3. 1-5) 
(EQUA3, 
STATEF) 

3, 1. 2 Nonlinear Aerodynamics Model 

The nonlinear aerodynamics model employs bivariate tables of lift and drag 
coefficients as functions of angle of attack and Mach number. These tables 
are fitted with bicubic spline functions to yield continuous first, second, and 
mixed partial derivatives during program execution. (BLYNE, BLICO, BLINE) 

3. 1. 3 Base Drag 

Base drag, D^, in the program is a function of altitude. It is an input uni- 
variant table. 


q = 1/2 P V 2 


3. 1. 4 Moment Balance Aerodynamics Model 

The moment balance aerodynamics model includes the effect of aerodynamic 
and propulsive moment balance trim forces on the overall applied loads. The 
distribution of trim between aerodynamic and propulsive is accounted for 
through a blend factor, j. The moment balance diagram for this model is 
given in Figure 3. 1-1. 

Asymmetric linear aerodynamics are employed with the addition of aero- 
dynamic moments. 


The uncorrected lift (untrimmed lift) is 


Y ( C L a + C L n ) 

u \ a O/ 


(3. 1-6) 
(AERO CO, 
BEROCO) 


(3. 1-7) 
(VT, UT) 


The equation for drag is 


: D = C D +k ( C L ) 
o \ u / 
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(3. 1-8) 
(AEROCO, 
BEROCO) 





In developing the equation for the total aerodynamic moment around the 
instantaneous vehicle center of gravity, the following assumptions are made. 

A. The base drag acts parallel to the vehicle axis and is centered at 
station Zg^. 

B. The axial contribution of aerodynamic trim force may be neglected. 
This assumption is very good for aerodynamically stable airframes 
with trim surfaces aft of the center of gravity. The assumption is 
poor for moderate to highly unstable or extremely stable airframes 
since in either case, induced drag changes due to trim deflection 
should then be accounted for. 

C. The aerodynamic lift coefficient slope, C T , should be fitted to 

1-1 a 

vehicle data without aerodynamic trim surface deflection 
(untrimmed). 

The untrimmed moment about the center of gravity is (approximately) 

-*CG = (L u cos a + D sin a) (X CG - X CGR> 


- (D cos a - sin a) (Z^^. - Z^^r) 


( C Ot +C ltf Q ) 

\ o a / 


S Ref d Ref 


d b (z bd 


Z CG> 


(3. 1-12) 
(VT, UT) 


To distribute the trim force between the aerodynamic trim device (tail) and 
the engine thrust gimbal contribution, the blend factor, j, is used, j is an 
input function of dynamic pressure, q. 


j = j (q) 


(3. 1-13) 
(EQUA3, 
STATEF) 


The tail contribution then is 


( X T " X CG^ AL T _ ^CG 
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Ox 


(3. 1-14) 



or 


AL 


T 


j^cG 

< X T - X CG> 


(3. 1-15) 
(VT, UT) 


The corrected or total lift then becomes 


L 


L 

u 


+ al t 


(3. 1-16) 
(VT, UT) 


The engine thrust contribution is 


(X E ' X CG* Sin 6 E " (Z E " Z CG^ COS 6 E 


T = (1 - j)^ CG (3. 1-17) 


The equations above are employed in developing the governing equation set, K, 
for the in-plane control vector, w. The solution procedure is described in 
Section 9.2. 


3. 2 PROPULSION MODELS 

Four propulsion models are employed in PADS. The first is a simple 
rocket model using input vacuum thrust with optional throttling and comput- 
ing fuel flow. The second, for use with the SSSP sizing module, simulates 
two engines with different I 's, with optional throttling. The third model is 
an air breather simulation. The fourth model has dual parallel-burn engines. 


3. 2. 1 Simple Rocket Model 

The rocket vacuum thrust per engine, ^yAC’ may be in P ut as a tabu l ar func- 
tion of burn time or as a constant. The net thrust of the vehicle may then 
be calculated as: 


T 


F 

_ VAC 


' a exit 



x t mult 


(3.2-1) 

(EQUA3, 

STATEF) 
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where 


a exit 1S 


the nozzle exit area, 


p is ambient static pressure, 
r a 


T MULT 


is the number of engines or thrust multiplier. 


The vehicle rate of change of mass may then be calculated. 


m 



T 

VAC 


I 

sp 


MULT 

g r 


(3. 2-2) 
(ACCEL, 
APPLY) 


When acceleration limit throttling is employed, the net thrust, T, is computed 
using an appropriate governing equation (K) so that the total vehicle accelera- 
tion is bounded. The total vacuum thrust may then be calculated 


T VAC T + A EXIT P a T MULT 


(3. 2-3) 
(TH, FH) 


If no I loss table has been input, 
sp r 


the vehicle rate of change of mass is then 


m 


T 

VAC 


I 

S P 


(3. 2-4) 
(ACCEL, 
APPLY) 


If on the other hand, an I loss table 

sp 


%I = %I 
sp sp 


VAC 


) 


/rated / 


xlOO 


x . 01 


(3. 2-5) 
(IMPUL, 
IMPULS) 
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has been input, the effective I will be used and 

Sp eff 


Isp eff 


I x %ISP 
sp 


/^VAC\ 
\ ^rated J 


x . 01 


and 


T 

VAC 
I g 
sp eff r 


(3.2-6) 

(IMPUL, 

IMPULS) 


(3.2-7) 

(ACCEL, 

APPLY) 


where F ra ^ is either input as a constant or interpolated from the input thrust 
table at the initiation of the thrusting (initial burn time). 


3. 2. 2 Dual Engine Model 

The dual engine model used with SSSP sizing problems first calculates the 
effective I of the two engines with ^^AC ant ^ ^VAC ’ ^eir rate ^ vacuum 

thrust levels. 


T 


t vac t vac 1 + t vac 2 

I sp 1 I sp 2 ( T VAC 1 + T VAC 2 ) 
sp T VAC 1 I sp 2 + T VAC 2 I sp 1 


(3. 2-8) 


(3. 2-9) 
(IS PRAT) 


This model assumes that only the second engine may have time-varying thrust 
(table input) or may be throttled. The net effective I for the system will 
change if this is permitted. We define A as the net nozzle area and given 

C 

the net instantaneous thrust T, we may calculate the ratio of actual vacuum to 
rated vacuum thrust of the two engines, 
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T + A p 
e *a 

T + T 

VAC 1 VAC 2 


(3. 2-10) 
(ISPRAT) 


We also define 


K 1 


I 

sp 2 



(3. 2-11) 
(ISPRAT) 


We may then calculate the ratio of instantaneous 
Equation (3. 2-9). 



to that defined in 


I 1 
sp 

I 

sp 


R' 


K' T + T 

VACj + VAC 2 


( K ' + R' - 1)T VACi + R- T VACz 


(3. 2-12) 
(ISPRAT) 


The equation for effective I and its partial derivatives are programmed in 

subroutine ISPRAT. The program logic is designed to parallel the I loss 

sp 

table results by using the same type of dependencies; that is, the result of 

both calculations is the percent I as a function of percent rated thrust. 

^ sp r 


3. 2. 3 Air-Breather Propulsion 

The air-breather propulsion model employs bivariate tables of thrust, T, and 
specific fuel consumption, SFC, as functions of relative velocity, V, and 
altitude, h. The SFC is in the units of fuel per lb-thrust per hour. The rate 
of change of vehicle mass may be calculated as follows: 


(SFC) T 

g 3600 
& r 


(3. 2-13) 
(ACCEL, 
APPLY) 


3. 2. 4 Parallel Burn Propulsion Model, JPR0 = 3, JAER = 3 

The configuration of the parallel burn model is shown in Figure 3. 2-1. In the 

equations, all aerodynamic coefficients are functions of Mach Number, M. 
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Aerodynamic lift (untrimmed): 

L u = ( C L a * +C L 0 ) ,SRef 
Drag coefficient: 

C D = C D o ^ ( C L J 

Aerodynamic moment: 

^a = ( C ^ 0 “ + C ^ o ) q S Ref d Ref 
Center-of-g ravity dependency: 

Z CG* = Z CG ^ 


( 3 . 2 - 14 ) 


( 3 . 2 - 15 ) 


( 3 . 2 - 16 ) 


( 3 . 2 - 17 ) 
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where W is vehicle total weight 


X CG X CG (W) 


(3. 2-18) 


For this model, 


base drag is centered at station Z__. 

dD 


The total untrimmed moment, «^ GG , about the instantaneous center of gravity 
depends on which engine is being gimballed or whether both engines are 
gimballed. Assume the (3 engine, having thrust T^, is not gimballed but has 
fixed engine deflection, 6^ . Then 

E P 

^CG = ( L u COS 0 + D sin 0 ) ( X CG " X CGr) 

- (d cos o' - L u sin a) ( Z CG - Z CGR ) + ^ a + D B ( Z BD “ Z cg) 

_ ( Z E “ Z Cg) T P COS 6 Ep + ( X E " X Cg) T (3 sin 6 Ep 

(3. 2-19) 

If both engines are gimballed, the last two terms involving Tp are excluded. 


The tail contribution to balance the moment depends on the blend factor 
j = j (q), where q is the dynamic pressure ^ 


al t 


j 


CG 


( x t- x cg) 


(3. 2-20) 


For one-engine fixed, the gimballable engine deflection for engine y required 
to balance the remainder of the moment is: 



- X 


CGy 


sin o 


E. 




cos 



d-jJ^CG 


(3. 2-21) 


If both engines are gimballed, we have 


[ T (3 ( X Ep " X Cg) + T, Y ( X E y " X Cg)] Sin 6 E 
‘ [ T (3 ( Z Ep " Z Cg) + T Y ( Z E y " Z Cg)] COS 6 
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E 


(l-j)^ CG (3.2-22) 



This completes the governing equation for engine-deflection solution. It 
remains to describe the applied-load terms in the equation of motion. The 
general form is 


where f (y) are the applied-load independent terms, G is a diagonal matrix of 
applied-load independent multiplying factors, and a is the acceleration vector 

having elements 


a 



( 3 . 2 - 24 ) 


ana 

diag [G] = ^1,‘y’ v cosy’ g r j 

The equations for the first three elements of a are 


( 3 . 2 - 25 ) 


T ^ cos 


V _ 


(« +6 El) 


+ T„ cos /a+ 6 


(‘ 


\ _ D - D„ cos a 

■) * 


( 3 . 2 - 26 ) 


m 



( 3 . 2 - 27 ) 


then 


v A 

a 1 = a cos <f> 

X . 

= a sin <p 


( 3 . 2 - 28 ) 
( 3 . 2 - 29 ) 
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The equation for 


a m depends on Tj and T^. (Note: may not be throttled) 


a 


m 




where 

Tj and are functions of t 


(3. 2-30) 


and 

ISP^ and ISP^ are constants. 

3. 3 ATMOSPHERIC MODELS 

Three atmospheric models are available in PADS: 

A. 1962 standard 

B. 1963 Patrick AFB 

C. Vacuum. 

3. 3. 1 1962 Standard Atmosphere 

The 1962 Standard Atmosphere model is a coded version of the analytic 
representation of Reference 2 up to 195 km. Above 195 km, a 50-point 
weighted least squares polynomial extends the range of data to beyond the 
altitudes of interest. This model is contained in subroutines ANLATM in 
the steepest descent portion of the trajectory module and in ANL62S in the 
quasi-linearization portion. 

3. 3. 2 1963 Patrick AFB Model 

The 1963 Patrick AFB model is described in Reference 3. It employs 14th- 
order polynomials for atmospheric properties up to 400, 000 ft. 
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Section 4 

AERODYNAMIC HEATING MODEL 


The aerodynamic heating model employed in PADS calculates the stagnation 
heating rate on a spherical nose cap. The equation is 


Q Q MULT 



(4-1) 

(DER3A, NLDRV) 


where the nose radius is assumed to be 1 and is the base density of the 

atmosphere model. Q MULT is an in P ut fla g havin g a value of either L or 
zero depending on whether the heat load should be calculated in that portion 

of the trajectory. ^mulT necessar ^Y an arc- bependent flag. 
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Section 5 

MISSION CAPABILITIES 


The mission capabilities of the trajectory program are described in this 
section. These mission capabilities are, in a mathematical sense, the 
boundary conditions on the multi-point boundary value problem. By way of 
introduction to these boundary conditions, some definitions will be useful 
in later discussions of the calculus -of- variations formulation (see 
Section 17). 


There are two types of boundary conditions, initial conditions and targets. 
Initial conditions apply to the beginning of an arc whereas targets occur only 
at the end of an arc. Several types of initial conditions can occur in a 
problem. These are 

A. Fixed initial condition 


y^ = (known value) 

where y^ is element of the state vector. 

B. Continuous initial condition 

^i |t+ ~ ^i |t- 

where | t denotes the arc end point 

C. Known or computable mass discontinuity 


(5.0-1) 


(5.0-2) 


m | T+ = m | T — Am (5.0-3) 

D. Mass distribution 

Let superscript 1 denote branch one and superscript 2 denote 
branch two. Then 

(5. 0-4) 


m 


T + 


+ m 
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The target conditions that can be met are described in the remainder of this 
section. These targets are divided into ascent, entry, or auxiliary boundary 
conditions. In addition to those described, it should be noted that the rela- 
tive state vector is also a set of candidate targets. It should also be noted 
that all targets described are functions of the relative state and time only. 

5. 1 ASCENT BOUNDARY CONDITIONS 

The orbital parameters are illustrated in Figure 5-1, programmed in PDBC 
and PDBCQL, and listed below. 



PLANAR ELLIPTICAL 


CR155-I 

P r SEMI-LATUS RECTUM 
R p PERIGEE 

R APOGEE 
a 

^ = R P + R a SEMI-MAJOR AXIS 
2 

R R 

e - a' p ECCENTRICITY 


H 

E 




= 2 GM/a, TRUE ANOMALY = TAN 
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A. Inertial velocity 



F. Orbital eccentricity 
Let 



G. Orbital inclination 

i = cos 1 (cos p sinijjj) 
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H. Argument of perigee 


3 = sin 

P 




I. Longitude of the ascending node 

= p 


/sin 4> -j- sin p\ 
I \ sin 1 f 


J. Semi-major axis 


R GM 


s 2 GM - RVj 


K. Apogee radius 

R = a (1 + e) 
a s' 


L. Perigee radius 


R = a ( 1 - e) 
p s 


M. True anomaly 


= tan 



The following three parameters are applicable to asymptotic injection, 
parameters are illustrated in Figure 5-2. 



B. 


Y = 


p e 
r r 

~ 2 
yjl-e 


The 


5-4 



CR155-I 


0 ASYMPTOTE 
Y NORMAL DISTANCE 
X HYPERBOLA AXIS DISTANCE 

X 


Figure 5-2- Asymptotic Parameters 


C. @ = cos * ^ ) 


Two additional orbital parameters are energy and momentum. 

A. Energy 

e = Lgm 
a 

s 

B. Momentum 

H = R Vj cos 7^ 


5. 2 ATMOSPHERE ENTRY MISSIONS 

Figure 5-3 illustrates the how total range, downrange, and cross-range are 
defined. 


CR155-I 


= REFERENCE LONGITUDE 
P f = REFERENCE LATITUDE 
j Jr = REFERENCE AZIMUTH 
U = VEHICLE LONGITUDE 
p = VEHICLE LATITUDE 
S D = DOWNRANGE 
S c = CROSS-RANGE 
S J = TOTAL RANGE 




Figure 5-3. Range Targets 






5. 3 AUXILIARY BOUNDARY CONDITIONS 

Additional boundary conditions in PADS for various purposes are given 
below. 

A. Dynamic pressure 



B. Heating rate 

® = °MULT 17600 ^( 26000 ) 

C. Reynolds number (unit) 

R =£ 
ey v 

The payload boundary condition equations are described in Sections 6 and 13. 
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Section 6 
PHASE I SIZING 


There are six sizing options in the Phase 1 sizing. These options are: 

1. Maximize payload or burnout weight with fixed initial weight 

2. Minimize lift-off weight with a fixed payload 

3. Minimize lift-off weight with a fixed payload and second stage 

4. Minimize lift-off weight with a fixed payload and first stage 

5. Maximize payload or burnout weight with a fixed (T/W)^ q and 

thrust 

6. Maximize payload or burnout weight with a fixed (T/W)^ „ and 
initial weight 


Each option is solved in its own subroutine and discussed below. 


6. 1 OPTION 1 

The equations for sizing option 1 are based on the initial weight and the mass 
ratios that come from the trajectory program. The booster propellant weight 
is given by 


W 

P B 



O. 



( 6 . 1 - 1 ) 


and the burnout weight by 


= W T _ - W 

BO L . O . p 


B 


( 6 . 1 - 2 ) 


The orbiter inert weight is determined from either a tabular input of the 
booster stage weight or from the following expression for this parameter 


W 

e 


B 


+ a 


w 


+ a, w 


1/3 


B 


B 


+ a_ w 
3 p 


2/3 


B 


(6. 1-3) 
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or 



f (W 


B 


Then, 


W 

o 


W 


BO 


- W 

e 


B 


The rest of the stage parameters are known. 


6. 2 OPTION 2 

The second option requires the sizing to be done from the top down, 
orbiter mass ratio is determined from 


= exp 


4V TOT ' e r Is PB <n 0*B)' 


g I 
e r sp. 


and the propellant weight is given by 



where 


W 

e 


o 


b + b- 


1/3 

W + b 0 W ' + b. 

p 2 p 

r o o 



or is input via tabular data. 

The orbiter initial weight is given by 


W = W + W + PL 
op e 

^o o 


(6. 1-4) 


The 


( 6 . 2 - 1 ) 


( 6 . 2 - 2 ) 


(6.2-3) 


(6.2-4) 
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The booster weights are determined from 




+ W 


>) 


W BO = W o + 


W 


B 


(6. 2-5) 


(6. 2-6) 


and W P is determined from Equation (6. 1-3) or tabular input. The 1 
B 

weight is given by 


W 


L. 


O. 


w 

B BO 


(6. 2-7) 


6.3 OPTION 3 

The option 3 sizing starts with the determination of the orbiter gross weight 

w = W + w + PL (6.3-1 

op e 

r o o 

These quantities are input. 

The booster mass ratio is given by 



The propellant weight is given by 


W 

PB 


C'B- 1 ) ( W o + W e B ) 


(6.3-2) 


(6.3-3) 
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where W e g is given by Equation (6. 1-3). The booster burnout weight is 
determined from 


W 


BO 


W 

o 


+ W 

«B 


(6.3-4) 


and the booster lift-off weight by 


W 


L. O. 


- W 


B BO 


(6.3-5) 


6.4 OPTION 4 

This sizing option requires the minimization of the lift-off weight. The 
booster lift-off weight is given by 


W 


L. O. 



(6.4-1) 


and the burnout weight by 


W 


W 


L. O. 


BO 


H- 


B 


(6.4-2) 


The orbiter initial weight is given by 


W 

o 


W 


BO 



B 


(6.4-3) 


where W 6 g is determined from Equation (6, 1-3). 

An iteration is required to determine the orbiter size to complete the mis- 
sion with the fixed payload. The following equations are solved iteratively 
until the payload error, € , is within tolerable limits. 


M- 


o 


EXP 


AV 


TOT 


- g I 
& r sp 

g I 
°r sp 
*o 



(6.4-4) 
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(6.4-5) 



W = W - 
P o 

l o 



W e = f /W p jor (6. 2-2) 


PL = W r - W 
f e 

o o 


PT 

-^FIXED 


PL 


The following test is made 


(6.4-6) 

(6.4-7) 

(6.4-8) 

(6.4-9) 


if PL < jxED* . P B ” ® go ^o (6.4-1) 

PL = PL fixed ; EXIT (6.4-10) 

PL > PL fix£d ; P b = P B + « go to (6.4-2) 

6. 5 OPTION 5 

This option takes advantage of the equations derived in Section 6. 1. The 
lift-off weight is determined from 


W 


LO 


N ( T VAC - A e Pa) 
(T/W) 


(6. 5-1) 


and Equations (6. 1-1) through (6. 1-4) are solved. 
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6. 6 OPTION 6 

Option 6 also takes advantage of Equations (6. 1-1) through (6. 1-4). The 
trajectory thrust is modified by 

T VAC =(#)Gtt) +A ePa f 6 - 6 - 1 ' 

6.7 PAYLOAD BOUNDARY CONDITIONS 

The payload boundary condition is available in the trajectory module of 
PADS for use with Phase I sizing problems. This boundary condition is 
employed for optimal staging problems (rubber stage). The equation for 
payload boundary condition is the same as Equation (6.4-8). The implica- 
tions of rubber- stage optimal staging are described in Section 13. 
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Section 7 

INTERFACE WITH SPACE SHUTTLE SYNTHESIS PROGRAM 


As was mentioned in Section 1 of this volume, the manned reusable launch 
vehicle synthesis is performed by the SSSP module in PADS. The synthesis 
modeling and equations remain essentially the same as documented in Ref- 
erence 1. There are, however, some minor model changes and program- 
ming changes. The model changes will be discussed below. The 
programming changes, which are mostly related to data communication, 
are listed herein and discussed in Volume II of this report. 

7. 1 THRUST SIMULATION CHANGE 

The thrust simulation available in the original SSSP program is completely 
dependent on its trajectory module (GTSM). In PADS, the thrust simulation 
is likewise related to the TABTOP trajectory module. This thrust simulation 
is described in Section 3. 2 of this report. The information on thrust and fuel 
flow values must be transferred to the SSSP program in order to calculate 
fuel weights. This calculation of fuel weights and related quantities is per- 
formed in a new subroutine called THRUST. 

7. 2 ENVIRONMENTAL EFFECTS ON DESIGN 

A new thermal protection system weight estimation equation has been added 
to the SSSP module. It has the form 

W TPS = (°l) 1 1 t L ) 2 C TP 1 + C TP 2 ("7.2-1) 

where C^Pj and C'pP 2 are i n P ut coefficients and r ^ and are input 
exponents. 0 L is the effective heat load (Btu/ft 2 ) on the entry portion of 
the trajectory. The effective heat load, 0 L> is computed by considering 
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the integral of Q 


only when Q 


is above an input threshold value, 0-p^ 



( 7 . 2 - 2 ) 


The quantity t^ in Equation (7.2-1) is 

‘L* t 2-*l ,7 - 2 - 3) 

meaning the duration when O is greater than 

A familiar configuration of coefficients for TPS design is 
1/8 
3/8 

(7.2-4) 

1 . 

0 . 


'TP. 


TP. 


7.3 HUNTING PROCEDURE 

A parameter hunting procedure for solving bounded optimization of up to 
10 design parameters used in SSSP design synthesis is available in PADS. 
The technique employed is called Powell's method, originally published in 
Reference 4. Bounding of free parameters is accomplished through the 
Box transformation of Reference 5. 


Briefly described, the method required no gradients and employs a conjugate 
direction quadratic ray search to find a minimum. The equations and pro- 
gramming logic for subroutine POWELL are presented in detail in Volume II. 
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7.4 MISCELLANEOUS CHANGES 

Following is a list of key programming changes in SSSP: 

A. Input data communication 

B. Merging of fly-back range calculations 

C. Traj ectory program communication 

D. Output format 

E. Sizing options for solids. 
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Section 8 

AUXILIARY PRINT COMPUTATIONS 


This section will detail the auxiliary computations required at each print 
point. They consist of the instantaneous impact points, inertial Cartesian 
coordinates, Euler angles, steering angles, and some orbit parameters. 


8. 1 ADDITIONAL ORBIT PARAMETERS 

The majority of the orbit parameters printed at each time point are calcu- 
lated as described in Section 5. Three parameters printed are not calculated 
in that section. They are the apogee and perigee velocities and the orbit 
period. They are determined from the following equations. 



( 8 . 1 - 1 ) 


( 8 . 1 - 2 ) 


(8. 1-3) 


8.2 INSTANTANEOUS IMPACT POINT 

The instantaneous impact point (IIP) is that point on a spherical earth where 
the vehicle would impact if it continued on its current path. The solution 
assumes unpowered vacuum flight on a Keplerian orbit. 

Certain quantities are calculated and tests are made before the IIP can be 
calculated in subroutine CRASH. The first test is made to determine if the 
orbit will intersect the earth. If the perigee radius is greater than the earth 
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radius, the orbit does not intersect the earth and a message is printed. If 
the apogee radius is less than the earth radius, an error has been made 
and a message is printed. 


The current true anomaly, t, , is given by 


p r /R - i 

cos 4 = 


( 8 . 2 - 1 ) 


The true anomaly at the impact point is given by 


COS ^ IIP 



( 8 . 2 - 2 ) 


The eccentric anomaly at the current position and the impact point are given 
by 


cos E 


e + cos £, 

1 + e cos £, 


( 8 . 2 - 3 ) 


sin E 


J 


i - e 


sin t, 


1 + e cos t, 


( 8 . 2 - 4 ) 


E 



( 8 . 2 . 5 ) 


The impact velocity is given by 


V 


IIP 



( 8 . 2 - 6 ) 


The impact elevation angle is given by 


cos 


R Vj. cos Y 

Y IIP = e r v iip 


( 8 . 2 - 7 ) 
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The true anomaly to the impact point is 





(8 


The geocentric impact latitude is given by 

sin Pjjp = cos P cos 4*2 s i- n ^ + s ’- n P cos ^ 

The azimuth angle at the impact point is given by 

cos p sin 
Sin ^IIP = cos p n p 

( cos i sinp^p - sin p \ 

= ) cosp IIP 

sin t, / 


(8 


( 8 . 


( 8 . 


* 


HP 


tan 


/ sin *IIP \ 

\ cos ‘•'iip/ 


( 8 . 


The longitude increment from the burnout point to the impact point is 


sin p 


sin C sin 4<j];p 

COS Pjjp 


( 8 . 


cos s - sin 

P sin P np 

(8. 

COS p - cos p cos 

P IIP 

? = tan" 1 ( SiniI ) 


(8. 

\ cos p / 




. 2 - 8 ) 

. 2 - 9 ) 

2 - 10 ) 

2 - 11 ) 

2 - 12 ) 

2 - 13 ) 

2 - 14 ) 

2 - 15 ) 
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The eccentric anomaly at the impact point is given by 


1 - e 


sin £, 


sin E 


IIP 


IIP 


1 + e cos £, 


( 8 . 2 - 16 ) 


IIP 


e + cos 


cos E 


IIP 


IPP 1 + e cos t, 


IIP 


( 8 . 2 - 17 ) 


E IIP = tan 


-i / Sin E IIP 


cos E 


IIP 


( 8 . 2 - 18 ) 


The time to impact is based on the difference in the time since perigee as 
determined at the burnout point and at the impact point. These times are 
given by Kepler's equation as follows: 


T q = X (E - e sin E) 


( 8 . 2 - 19 ) 


and 


T IIP = X (E IIP " ® Sin E IIP } 


( 8 . 2 - 20 ) 


where 



The time to the impact point from burnout is thus given by 



( 8 . 2 - 21 ) 


( 8 . 2 - 22 ) 
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t 



and 


r = r TTX , - T + r if T < o. (8.2-23) 

IIP o 

The impact longitude is given by 

~HP = (8.2-24) 

The impact time is given by 

T = T + T (8. 2-25) 

1 IIP o v 

The range to the impact point is 

S = E r T (8.2-26) 

The range from the reference latitude and longitude point is given by 

cos X = cos p R cos p IIp cos (p R - M-jpp) + sin P R sin P n p 2-27) 

S np = E r cos ’ 1 (cos X ) (8. 2-28) 

The azimuth angle from the reference point to the impact point is given by 



(8. 2-29) 
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However, if 


(p R - n IIP ) = 0 and 


P HP < P R : A *IIP = w 


P IIP - P R : A "lIP ° 


(8. 2-30) 


or, if 


(h*R - H np ) = v and 


P IIP < " P R : A>,< IIP V 


P > - P • A * =0 

H IIP - P R- HP 


These quantities define the attitude and position of the vehicle at the impact 
point. 

8. 3 INERTIAL CARTESIAN COORDINATES 

The inertial Cartesian coordinates are determined from the inertial longi- 
tude, the current latitude, and the radius to the vehicle. The Cartesian frame 
has the following as its primary planes; the equatorial plane, the plane through 
the launch site and the north pole, and a plane perpendicular to these planes. 
The position transformation is given by 


Xj = R cos p cos p.j. 


(8.3-1) 


Yj = R cos p sin p. (8. 3-2) 

Zj = R sin p (8. 3-3) 


The velocity transformation is obtained by differentiation of the above equa- 
tions to give 

• • • • 

Xj. - R cos p cos fij - R p sin p cos p^ - Rp^. cos p sin p^ (8. 3-4) 
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(8. 3-5) 


Yj = R cos p sin p ^ - Rp sin p sin p^. + Rp^ cos p cos p^. 

• • • 

Zj = R sin p + R p cos p 

where 
R 
P 

• V cos Y sin di 

p _ X. 

R cos p 


= V sin V 
V v 

= — COS Y COS 


(8. 3-6) 


(8. 3-7) 
(8. 3-8) 

(8. 3-9) 


8.4 RELATIVE EULER ANGLES 

The relative Euler angles of the vehicle axes are obtained from the vehicle 
state by establishing a Cartesian coordinate system at the vehicle center of 
gravity and performing the following rotations in order: 

1. Rotate about the z axis through the azimuth angle kjj . 

2. Rotate about the y axis through the flight path angle Y . 

3. Rotate about x through the bank angle <{> . 

4. Rotate about y through the angle -of -attack, a . 


The matrix product of these rotations results in the matrix of direction 
cosines between the vehicle axes and the relative Cartesian set. 


D 


cos a 0 sin a 


( Direction \ 
Cosines 


1 0 
0 cos a 




(8.4-1) 


Given the direction cosines, the Euler angles are calculated in subroutine 
DCTOE for both yaw-pitch- roll and pitch-yaw-roll sequences of rotation. 
The equations for these angles are YAW, PITCH, ROLL 


9 PITCH = tan 


-1/ D (1, 3) 


V >/l - D ( L 


3)‘ 


(8.4-2) 
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S YAW = tan 


0 


ROLL 


1 ( P ( 1 , 2 ) \ 

( 1 , d) 

-1 ( D (2, 3) \ 
\ D (3, 3) J 


(8.4-3) 

(8.4-4) 


PITCH, YAW, ROLL 


0 


PITCH 


tan 


1 / -D (1, 3) \ 
V D (1, 1) ) 


(8.4-5) 


0 YAW = tan 


-1 / -D(l, 2) 


0 


ROLL 


tan 


1 - D (1, 2) 


-1 / -D (3, 2) \ 

\ D (2, 2) j 


(8.4-6) 


(8.4-7) 


where D (I, J) is a location in the direction cosine matrix given in 
Equation (8. 4-1). 

8.5 INERTIAL EULER ANGLES 

The inertial Euler angles are obtained by a right multiplication of the relative 
direction cosine matrix [(Equation (8. 4- 1 )] by the transformation 



(8. 5-1) 


to obtain the matrix of inertial direction cosines 


/ direction\ _ / direction \ , D ,, 

\ cosines ) \ cosines / 


(8. 5-2) 


R 


The same equations used to determine the relative Euler angles (subroutine 
DCTOE) are used to determine the inertial Euler angles from 
Equation (8. 5-2). 
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8, 6 STEERING ATTITUDE ANGLES 

Two steering attitude angles are calculated for output. The first angle is the 
elevation angle from the local horizontal to the vehicle centerline in the 
plane determined by the local vertical and the velocity vector. The second 
angle is the azimuth angle measured from north to the vehicle centerline 
projection on the local horizontal plane. The equations for these angles 
are: 


0p = Y + sin ^ (sin a cos <j>) 


0 - 1 
l = 4* + sin (sin a sin 4> ) 


( 8 . 6 - 1 ) 
( 8 . 6 - 2 ) 
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Section 9 

CONTROL VECTOR SOLUTION - NON-OPTIMA L CONTROL 


To integrate the equations of motion it is necessary to know the control 
vector 


U 


(T, 6 e , 


a, 4> ) 


(9.0-1) 


Because we are concerned with endo- as well as exo-atmospheric trajec- 
tories, this can be a formidable problem, even if we overlook the added 

complication of determining the optimal steering angles, a and <}> . (See 

Sections 11, 12, 16, and 17 for discussion of this problem. ) 

For example, suppose we want the vehicle to perform a "gravity turn' while 
it balances its aerodynamic moment and, at the same time, conforms to a 
total acceleration limit (not an unlikely circumstance). Under these condi- 
tions, none of the in-plane control quantities, T, & E> or a, can be explicitly 
evaluated; each is dependent upon the other two. In other words, we are 
confronted with a situation that requires the solution of a system of nonlinear 
simultaneous algebraic equations. Were this the only example in which this 
predicament arose, it would be of little burden, but, unfortunately, it is 
only one of many. 

The difficulty of the problem is compounded by the fact that the governing 
conditions may change at any time. For example, the thrust may have to be 
instantaneously throttled in order to satisfy a total acceleration limit, or 
the angle of attack may have to switch from one mode of flight to another. 

Hence, to determine the vector U, the above difficulties compel us to adopt 
the rather oblique and formalistic approach described below. In doing so, we 
reap at least one important supplementary benefit - later program modifica- 
tions will be greatly facilitated by the formalism we introduce now. 
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9.1 NON-OPTIMAL BANK ANGLE 

Presently, there are only two non-optimal conditions for the bank angle, 4> „ 
One of them is that 4> = 0. The other is the vertical rise/pitchover mode. 

In both cases, an explicit equation for <f> can be obtained. In the former 
case it is simply 

4> = 0 ' (9. 1-1) 

and, in the latter, it is 

4> = tan" 1 (k^/k^) (9. 1-2) 


where 


k ^ = V cos Y sin p 


V 

R cosp 


cosY sin 


i|i + 2 co ^ 


+ co cos p (R <0 sin p sin c)j - 2V cos c|j sin Y ) 


(9. 1-3) 


and 


k Y = (“IT " g) cos Y - v Y * 

+ co cos p ^2V sin v)j + R w (cos p cos Y + sin p cos c|> sin Y )J 

(9. 1-4) 

Asa result, we can exclude this component of the U vector from the formal- 
istic treatment and restrict our attention to the problem of determining the 
in-plane control vector 

w = (T, 6 E , a) T (9. i_5) 
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9. 2 IN- PLANE CONTROL DETERMINATION WHEN a IS 
NON -OPTIMAL 

Suppose the vehicle is at some arbitrary point of a non-optimal portion of its 
flight. Let 


K = (K r K 2 , K 3 ) T (9.2-1) 

be the system of governing equations for the choice of w at this particular 
point. In general, K will be a function of the state and the in-plane control; 

1 « G • y 


K = K (y, w) (9. 2-2) 


The in-plane control vector (value of w) will be determined by solving: 


K (y, w) = 0 (9. 2-3) 


To do so, we make use of the well-known Newton-Raphson iteration, which 
runs as follows: 


Starting with some initial guess for w, e. g. , w Q> compute an increment 
Aw from the equation 


(9.2-4) 

(ALGCON, 

BLGCON) 


where K (w ) is the matrix of explicit partials of K with respect to w 
w o 

evaluated at w = w ; i. e. 


K (w ).. = K. (w ), i - 1, 2, 3, ; j - 1, 2, 3 (9. 2-5) 

w o 11 i o 

J w. 

J 


Aw = -K” 1 (w ) K (w ) 
w o o 
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If 

(9. 2-6) 

(see Note) (ALGCON, 

BLGCON) 

we are done; otherwise, increment w q by Aw and repeat the process until 
Equation (9.2-6) is satisfied. 

A rigorous discussion of the necessary conditions for the convergence of this 
iteration would be inappropriate here. Suffice it to say that its success is 
largely dependent upon the guess w q and the behavior of K^. 

Presently, PADS' policy in choosing w q is the following: 

At the first point of the trajectory, 

w (t ) = (10 6 , 0, 0) T (9.2-7) 

o o 

at all corner points, 

w Q (t^) = w(t ) (9. 2-8) 

and at the interior points of a subarc 

w (t) = w (t - At) (9. 2-9) 

o 

where At is the integration step size and w is the converged value of the 
in-plane control vector. So far, this policy has been adequate, but the need 
for a more elaborate one can always arise. 

Note: « = 10 ^ CDC 6500 version 

_ 7 

10 Univac 1108 version 


3 
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9. 3 PARTIAL DERIVATIVES OF THE NON-OPTIMAL CONTROL 
VECTOR 

Both the steepest descent and QL modules of PADS require the evaluation of 
the total first partials of the in-plane control quantities at each point of the 
trajectory. Steepest descent needs them on its backward integration in order 
to construct the adjoint coefficients (see Section 12). The QL module needs 
them to evaluate the Euler-Lagrange equations (see Section 16). 

For non- optimal portions of the flight, it follows from the implicit function 
theorem that these partials are given by 


D 0 w 

y 



w 


K 

y 


(9. 3-1) 


where D is the total partial derivative operator 


D = 1 

y 1 

( 9 . 

9 

•• 9 ) 

(9. 3-2) 

1 9 V 

ay 2 ' ' ‘ 

a v 



K is the matrix of explicit partial derivatives of K with respect to w, and K 
w 

is the matrix of explicit partial derivatives of K with respect to y. Of cours 

both K and K are evaluated at the converged value of w. 
w y 


Because the QL module of PADS employs the method of quasi-linearization 
to solve the multipoint boundary value problem that arises from the calculus 
of variations, it also requires the total second partial derivatives of w with 
respect to y. Differentiating Equation (9. 3-1) with respect to an arbitrary 
state variable y. yields 



(Dy * w) 



w 


Kyy. + Y.+/K + X. \ D * w 

yy i i ^ wy. i ) Y 

i = 1, 


n (9.3-3) 


9-5 


n> ^ 



where the j-th column of the matrix Y^ is given by 


Y.^ = K D • w, j = 1, n (9. 3-4) 

i wy. Y i 


and the k-th column of the matrix X i is given by 


X. (k) = K Dy. • w, k = 1, 2, 3 (9.3-5) 

1 ww, 1 
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Section 10 
CONTROL LAWS 

In Section 2.4, the concept of the control vector U and its subset, the 
in-plane control vector, 

w T = (T, &£, a) (10.0-1) 

was introduced. From prior descriptions of the simulations in PADS and 
discussions in Section 9 , it is apparent that w can be calculated explicitly 
only in special circumstances. In this section, the equations corresponding 
to in the algebraic constraint vector introduced in Section 9 will be 
presented. The first sections describe the non-optimal control capabilities 
in PADS. The later sections describe the calculations necessary for 
in-flight bounded control and state function constraints. 

10.1 VERTICAL RISE AND PITCHOVER 

During the vertical rise or pitchover control mode, the bank angle and angle 
of attack are fully determined. The equation for bank angle is given in Equa- 
tion (9. 1-2). The equation for involves k y defined in Equation (9. 1-4) is, 

= (T sin (a - 6-g) + L - D^ sin a) cos <j> + m k y (10. 1-1) 

This equation plus associated partial derivatives may be found in subroutine 
BL4 (steepest descent) and A L4 (Q. L. ). 

10.2 CONSTANT ANGLE OF ATTACK 

The constant angle-of-attack control mode is used for a number of different 
situations in PADS. This simplest case is when the zero angle-of-attack 
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control mode is chosen. The form of constraint function, K^, in all such 
cases is, 

K = a - C =0 (10. 2-1) 

3 a 

The zero a control mode has 


C = 0 (10.2-2) 

a 

The most significant use of this control mode is in the steepest descent 
control computation. As shown in Equation (12. 1-23), the control correction 
6a is computed using the steepest descent equations. This increment is 
added to the nominal control 


a = a OLD + 6a (10 “ 2_3) 

This new value a is inserted into Equation (5. 2-1) by letting 

C a = a (10.2-4) 


This control constraint is likewise used under any circumstance when the 
maximum angle -of-attack limit is reached. Suppose a* is the calculated 
angle of attack either in the steepest descent or quasi-linearization program 
and 


a* > a (a > 0) 

max max 


(10. 2-5) 


Then 


C = 
a 


a 

max 


(sign (a*)) 


( 10 . 2 - 6 ) 


10-2 



10.3 GRAVITY TURN CONTROL 

The gravity turn control establishes a balance of aerodynamic and propulsive 
forces so that there is no net force normal to the flight path. Referring to 
Equations (2. 3-2) and (2. 3-3), this will occur when 

X 

T sin (a - & E ) + L - sin a = 0, or a =0 (10. 3-1) 

The form of the constraint is the same. 

I 

X 

= T sin (a - & E ) + L - sin a = 0, or K 3 = a =0 (10. 3-2) 

10.4 MAXIMUM LIFT CONTROL BOUNDARY 

The maximum lift force magnitude may be instantaneously constrained by 
employing this control mode. From Equation (3. 1-6), the uncorrected lift 
is 


L 

u 


( Cl o + C L °) q 6Ref 


for aerodynamic option 1 or 3 

(10.4-1) 


and 


L u ' C L q 6 REF 


for aerodynamic option 2. (10.4-2) 


The desired maximum lift 


C = L sign (L ) 
u T max 6 u' 

.Li 


(10.4-3) 


Hence, the form of the K 3 function is 


K, = L - C =0 
3 u u L 


(10.4-4) 
(A L3, BL3) 
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10. 5 UNPOWERED TOTAL ACCELERATION LIMIT 

When the unpowered total acceleration limit is applied, the angle of attack 
may be solved for, and the bank angle may still be optimized. The equation 
for is: 

K 3 = U + D 2 - (g max W) 2 = 0 or K 3 = (a*) 2 + (a V ) 2 - (g max W) 2 

(10. 5-1) 

This equation and it associated partials are programmed in . subroutine s AL5 
and BL5. 


10.6 STATE INEQUALITY CONTROL MODES 

When state instantaneous inequality functions are imposed on the trajectory, 
first of all, a corner point must be included at which the bounding function 
reaches its maximum value. After reaching its maximum value, the vehicle 
will fly the boundary until the optimized control will tend to tangentially fly 
the vehicle off the boundary. While the vehicle is on the boundary it is 
actually flying, a non- optimal control law which must be included in the 
set. 


There are three state inequality bounding control laws in PADS. These are 
listed below with their corresponding inequality limit. 


Control Law 

q = o 

Q =0 

R =0 

ey 

For the dynamic pressure rate, q, 


Corresponding Inequality 

Maximum dynamic pressure limit 
Maximum heating rate 
Maximum Reynolds number 

control law, the equation is 


. or a 

K_ = V R -rrr- + 2p Y = 0 
3 3h r a 


( 10 . 6 - 1 ) 
(A L7, BL7) 
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For the heating rate limit, Q, control law, the equation is 


. dp 

1 1 5 R — — 

K 3 - IV ' l.T ^ = » 


( 10 . 6 - 2 ) 
(AL8, BL8) 


For the Reynolds number rate, R , control law, the governing equation is 

ey 


dP, 


3 v 


K_ 


v 3h 9h ^a 


R + V = 0 


(10. 6-3) 
(AL9, BL9) 


where v is the viscosity. 
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Section 11 

POSING THE PROBLEM FOR STEEPEST DESCENT 

The steepest descent method as applied to solving for the control vector U 

will be described in the following steps. 

A. Derive the influence functions for steering angle, arc time duration, 
and initial state effects on arbitrary end boundary conditions. 

B. Using these influence functions, develop the equations for the closed 
loop steepest descent algorithm. 

C. Show how the closed-loop steepest descent algorithm generates 
steering angle and trajectory parameter corrections to first force 
the resulting trajectory to satisfy problem constraints and second 
drive a payoff quantity to its approximate optimal value. 

D. Show how the steepest descent solution may be transformed to 
approximate the exact Lagrange multipliers for use in starting the 
quasi-linearization module of PADS. 


11. 1 ADJOINTS AND INFLUENCE FUNCTIONS 
The equation 

y = f.(y, U., t. } 0 £t. ST. <U. 1- 

defines the dynamical system during arc i subject to an arc cut-off function 


fl. ( y, w., t. } = 0 ' 1 i . A " 

i l i i J 

which determines i\. The control vector, U, contains the optimizing steer- 
ing vector u (see Note). 


U. 

i 



(11. 1-3) 


Note: The solution for U on arcs with suboptimal control (K(U, y) = 0) is 

explained in Sections 9 and 10. 
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Then, the optimization problem resolves into determining a set of U histories 
(one per arc) and u . adjustable parameters that will optimize some function 
of state $ at the terminus of an arc while constraining additional functions, 

, at the termini of that or other arcs. 

If the trajectory during arc i as defined in Equation (11. 1-1) is perturbed, 
the perturbed solution can be related to the nominal result at any time t 


y' = + 6Y 


(11. 1-4) 


Discarding higher-order terms, this may be expressed in operator notation 
as a Taylor series which, after subtracting the nominal value, becomes 

a r. 

6y ' = D* • f. 6y + fiu. (11. 1-5) 

1 y i 7 3u. i 


where D 


y 


T, and 0f /3u^ 


are matrices of partial derivatives (See Note). 


The end conditions for the i th stage are given by the equation 



• ( 11 . 1 - 6 ) 


Equation (11. 1-5) gives the perturbation effect of the state and control on the 
derivative y'; however, we really need to know the effect of state and control 
changes on the end condition 


Note: 


Dy is a partial derivative operator which assumes the steering con- 
trol vector elements are constant. 
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We postulate- the desired form of this relation as 


6^. 

1 




+ E ( 6 u i} 


(11. 1-7) 


where the matrix of influence functions \ 1 relates the effect of changes in 
the state at time t on 6'!''^. The function E relates the effect of control 
changes up to time t on 6'3/^. For the sake of this discussion , assume final 
time is fixed and under this stipulation take time derivatives of Equa- 
tion (11. 1-7), resulting in 


6'k . 

i 


= 0 




6y + E 


( 11 . 1 - 8 ) 


Now, substitute Equation (11. 1-5) into (11. 1-8) 



(11. 1-9) 


or 





/ *.\ T 9f. 

6y + \ X / "au" 6u i + E 


0 ( 11 . 1 - 10 ) 


This equation must hold for any 6y and therefore the coefficient of 6y must 
equal zero 



( 11 . 1 - 11 ) 
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or in more familiar form, take the transpose to get the adjoint differential 
equations (see Note) 



( 11 - 1 - 12 ) 


Equation (11. 1-12) is used to solve for the influence coefficients subject to 
the boundary condition 




(t^ fixed) 


(11. 1-13) 


The last two terms in Equation (11. 1-10) may be used to solve for the control 
influence function 



(11. 1-14) 


If the equation is integrated backwards from to t, then 


E 


t 




(11. 1-15) 


There can be no influence of control changes at 


and the refore 


T. 


1 


0 


Note: On arcs where control is subject to the algebraic constraint, 

Equation (11. 1-12) becomes 


* . / . \ T * . 

\ = -ID* • f. - D • f. D • K (D • K)' 1 ) k 1 

\y i u iy u ) 
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Using the following definition for the impulse response function 


A 1 



( 11 . 


Equation (11. 1-15) becomes 


E 


t 




6u. dt 
l 


( 11 . 


If the final arc time, Tp is allowed to vary, the perturbation in the end 
dition may be described by combining results in Equation (11. 1-7) 


6 *. = 
l 




A 1 6u. dt + dT. 
i 


(H. 


Equation (11. 1-18) is a vector equation which contains an unknown final 
variation, dT^. However, the scalar cutoff function ^ defined in Equa- 
tion (11. 1-2) should determine the final time variation. Since Equa- 
tion (11. -1-18) is a general form, we may describe the variation of the 
scalar cutoff function in a like manner. 

6fl. = (\ l ) 6y 

Since £2 = 0 and 6fi. = 6u„ dT. may be solved for as follows: 
i l l i 


dT. = -r- 

/ £2A T 

6w. - \k / 6y 

(° S2, 

- / A 1 6u.dt 

1 

i 

t = 0 T 

i 

m 

i 


r 

T. 


A 


n. 

i 


6u.dt 

l 


+ £2df. 


( 11 . 


1-16) 


1-17) 


con- 


1-18) 


time 


1-19) 


1 - 20 ) 
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Defining 


V. £ 2 . 

X 1 = X 1 - -A X 1 (11. 1-21) 

£ 2 . 

l 

and substituting Equation (11. 1-20) into Equation (11. 1-18) yields the relation 


6 SE' 


i 




A 6udt + — — 

& 601 


( 11 . 1 - 22 ) 


The term x'^i satisfies the adjoint differential equations if its boundary 
condition is 


*£2. 



1 

a £2. 

i 

A. 

t. “ a y 

T. 

ft. 

8y 


i 

1 

1 



(11. 1-23) 


This equation is programmed in subroutine ADIC operating in conjunction 
with subroutine PDBC which supplies the matrix for nonlinear 

functions. 


11.2 MULTI- ARC PROBLEMS 

Section 11. 1 has developed the influence functions for a single arc problem. 
In this section the formulation is extended to multi- arc problems with only 
terminal constraints (Reference 6). 


Thus, in an i stage problem 

dtf = 6¥. 

i 



( 11 . 2 - 1 ) 


means that all constraint perturbations are determined at the end of the i*h 
arc. Now examine the effect of control and parameter changes in the last two 
arcs of the trajectory on the constraint perturbation. 
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= 


/ *.n. ,\ T r° >&.n. , r° 

= \ X. 1 1 / 6y | o + / A 1 1 6u^ j dt + / A 6Uj, dt 



T. 

1- 

l 



T. 

1 

l X 1 11 

) T y 


4- fw.i 



S-i 

Vi 

r * U u) . 

1 

n i 

T. 

1 



(11. 2-2) 


Following the form of Equation (11. 1-23), the boundary condition on the 
adjoints at the next to last arc corner point is: 




T 


i-1 


i i 


i-l 



(11. 2-3) 
(ADID3A) 


noting that 



(11. 2-4) 


for correspondence with (11. 1-23). 

Equation (11. 2-2) may be put in completely general terms by summing over 


the index j 

* r ° 

i-l 

f <S.O. ,,\ T 6io. 

d * = y [ 

A 1 ^ 6u.dt + / 

U 1 J+1 ) y , J 

j 

J *-• 

n. 

j=i t j 

j = l 





£ 2 . 


6oo 


T. 

1 


, * ( “)’ 


6y 


(11. 2-5) 
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which is the general form for control and parameter sensitivities with bound- 
ary conditions defined only at the trajectory terminus. 

Equation (11. 2-3) is programmed in subroutine ADID3A. Arc time sensitivi- 
ties are programmed in subroutine STAU. 


11. 3 MIXED BOUNDARY CONDITIONS 
A more general class of problems arise when constraints are not completely 
defined at the trajectory terminus; that is, d# 


The simplest example of this is the case where the boundary condition is an 
explicit function of arc times 


i 



(11. 3-1) 


We refer back to Equation (11. 2-4), which now becomes 



(11. 3-2) 


The equation for the particular adjoint discontinuity corresponding to Equa- 


tion (11. 2-3) 

x 11-1 

is 

't.n. 
x 1 1 


r m d$. 

Ax 11 ) y+a Vi 

a £2. , 

i - 1 

ay 


T i-i 


+ 

T i-1 

n i-i 



(11. 3-3) 
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For the case where 


I 6T j ' 


(11. 3-4) 


This corresponds to three different boundary conditions in PADS: 


A. Elapsed time 


- I 


B. Inertial longitude 


M-t = V- + 


C. Longitude of ascending node (see Section 5). 

Equation (11. 3-3) is programmed in ADID3A. 

11.4 INTERMEDIATE BOUNDARY CONDITIONS 

Suppose some components of the d* vector are completely defined at an arc 
corner point - say, arc s - prior to the i th arc. This implies that the ^ 
boundary condition vector can be partitioned 


¥ = (* s 


For the partitioned terms, Equation (11. 2-5) becomes 


£ f 6ujdt + 2 (.Vjllf y |1 


j = l T J 


J T. " 

J J 


/o 


* n. 

+ X S 1 6y 
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Likewise, the boundary conditions on the intermediate constraint adjoints are 


n 

V s s 


aT s ' 

_ is 

9D 

s 

ay 

T T 

Q 

o 

ay 


(11.4-3) 

(ADICB3) 


whereas the same arc corner the remaining adjoints are adjusted as in 
Equation (11. 2-3) or (11. 3-2). 



(11. 4-4) 
(ADICB3) 


These equations are programmed in subroutine ADICB3. 

11.5 BRANCHED TRAJECTORY BOUNDARY CONDITIONS 

Figure 11-1 illustrates the arrangement of a typical branched trajectory. 


t s defines the branch point and vector \& is partitioned into ^ and T k . 
Clearly all parameter and control changes daring segment I influence both 
and However perturbations during segment II can have no effect on 

and vice versa. 


The adjoint boundary conditions are: 
At (segment III) 




tf.fi. 
1 1 




9% 


*i 8 n i 


9y 


" b. 9y 


r . 
l 

T . 
1 

1 

T . 
1 

egment II) 



_ £ 

k 

*k 

dy 

T v 

r 

** * 


(11. 5-1) 
(ADICB3) 


(11. 5-2) 

(ADIC) 


Equation (11. 5-1) is solved and the adjoints then associated with tf^ are 
integrated backwards from 14 to t+. This integration is stopped and then 
Equation (12.5-2) is solved and the adjoints associated with tf k are integrated 
back to t+ at which point the ^ and tf k vectors and adjoint matrix are 
merged. Adjoint discontinuities (if appropriate) are calculated using Equa- 
tion (11. 3-3) or (11. 2-3) and the integration may proceed to initial time. 

The equations for the branched adjoint boundary conditions are contained in 
subroutine ADICB3. 
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Section 12 

CLOSED-LOOP STEEPEST DESCENT ALGORITHM 

In the prior section, the influence functions for control, initial conditions 
and staging parameters have been derived. In this section these influence 
functions will be incorporated into the closed-loop, steepest descent algo- 
rithm formulation. 


12. 1 USE OF INFLUENCE FUNCTIONS 

We repeat Equation (11. 2-5) making the following substitutions: 



( 12 . 1 - 1 ) 


Note that £2 = 1 for a timed arc cut-off; hence Equation (12. 1-1) becomes 



( 12 . 1 - 2 ) 

(STAU) 


Therefore, using the new notation 



A 1 &Ujdt + 


' 'If? 

y s viTj 

L t j 

j=i 


(12. 1-3) 
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We keep in mind that the vector may be partitioned for branching or 

intermediate constraints and now also stipulate that the last element in this 
vector be the performance index, d$, as shown below: 


d'J' 



d'k. d<$ 


(12-1-4) 


We further simplify the notation of Equation (12. 1-3) by combining free initial 
conditions and staging time sensitivity terms 


S 






6p 




where 



k is a scalar relation. 

y 

t 

o 


(12. 1-5) 


( 12 . 1 - 6 ) 


(12. 1-7) 
(SINIT) 


Equation (12. l-3)then becomes 


dty. 

i 



by 


o + 


2 


j=i 



'f.n. 

j\ 1 ^ 6ujdt + S 1 6p 


( 12 . 1 - 8 ) 


It is apparent that for a given finite order d^ vector there are an infinite 
number of solutions to Equation (12. 1-8); moreover, if particular sensitiv- 
ities are small, the associated contributions of parameter or control pertur- 
bations may be negligible relative to more sensitive contributors. To 
account for the discrepancies in sensitivities we assign a weight to each of 
the control and parameter sensitivities. 
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Let W be a diagonal control sensitivity weighting matrix. For this program, 
it has two elements: 



(12. 1-9) 


Let Y be a diagonal parameter sensitivity weighting matrix. 



(12. 1-10) 


To account for the multiplicity of solutions to Equation (12. 1-8), we will seek 
the one that forces perturbations along the direction of steepest descent to 
satisfy all the constraints. To do this, the control and parameter perturba- 
tions must be maximized. 


The metric of control and parameter perturbations given below is to be max- 
imized according to the development of Reference 7. 


(dP) 2 


I 





Y6p 


( 12 . 1 - 11 ) 


To the (dP)^ metric we adjoin the problem constraints through the Lagrange 
multipliers |i™. Using variational calculus, we will maximize the resulting 
functional given below. 
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S Sp 


(12. 1-12) 


Substituting Equation (12. 1-11) into Equation (12. 1. 12) and combining terms 
gives 

i 


■ I 


J = i 


J (bajwfaj - q T A 1 ^Sujjdt 



/ ~ \T 

“ 

ryi 



+ p 

d^ - y X J by 




to . 


T T ^i 

- 6 p Y 6p - |jl S 6p 


(12. 1-13) 


In order for J to be maximized (stationary), 6J must equal zero. We there- 
fore take the first variation. Note that W, Y and p are constants. 


6 J 



^6ujw6(6m) + 6 (du^wSm 


T 


H- 


A 1 J 6(6uj) jdt 


- 6p T Y 6(6p) - 6(6p) T Y 6p - p T S* l 6(6p). (12.1-14) 

T T 

The terms 6(6u )W6u and 6(6p )Y6p are equal to their respective trans- 
poses since both W and Y are diagonal matrices. 
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(s(6u T )W6u^ 


= 6u T W T 6(6u) 


But = W and therefore 


) T 

= 6u T W6(6u) 


Hence 


^6(6u T )W6u^ T = 6u T W6{6u) 


(12. 1-15) 


and likewise 


6(6p T )Y6p = 6p T Y6(6p) 


(12. 1-16) 


Substituting Equations (12. 1-15) and 12. 1-16 into Equation (12. 1-14) yields 


6 J 


■I 

j = l 


/•( 


26u T W6(6u.) - p T A 1 J f(6u.)jdt 
J J J 


T. 

J 


-^26p T Y + p T S^ 1 )t(cp) = 0 (12.1-17) 

The coefficients of 6(6u); and 6(6p) must equal zero for arbitrary 
perturbations. Thus 




1/2|J. T A 


*.S2. . 

1 %■' 


or 6t \i 


T ^i - 1 

l/2(i S Y or 6p = 



(12. 1-18) 


(12. 1-19) 
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The substitution of Equations (12. 1-18) and (12. 1-19) into Equation (12. 1-8) 
will allow the solution for 



- 1/2S 




( 12 . 1 - 20 ) 


or 


M- = 



tf: 


- S 



T 


-1 


X 





(12. 1-21) 


or, if the bracketed term is called the A matrix, 



(12. 1-22) 


Substituting Equation (12. 1-22) into Equations (12. 1-18) and (12. 1-19) gives 




(12. 1-23) 
(MTXI) 
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and 


) 12. 1-24) 

to / (MTXI) 

If the A matrix is inverted only at the initial point, in Equation (1 2. 1-23), 
the time history of impluse response functions will yield a history of control 
corrections which will under normal circumstances drive the d4f vector to 
zero. However, it is advantageous to consider that as the new trajectory is 
integrated, the difference in state between the nominal and new trajectory at 
each arc time point does affect the end constraints, d^. It is from this con- 
sideration that the closed-loop steepest descent algorithm arises (Refer- 
ence 7). In this algorithm, the time history of the control contribution to A 
is computed, and the parameter sensitivity contributions are included up to 
the time where the parameter is adjusted. An example of this is given below 
where the control correction 6u, at arc time t during arc k is shown 




6y (t} = y(t) - y(t) < 12 - !- 26 > 

This equation and the corresponding one for parameter corrections are com- 
puted in subroutine MTX3A. 

As described in Section 11.4 and 11.5, the d'T vector may be partitioned for 
intermediate constraints or branching. The A matrix is symmetric and of 
the same order as the number of problem boundary conditions. If during a 
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portion of the trajectory only part of the vector applies, the rows and 
columns of A associated with the active constraints are compressed into a 
new matrix, called B, before inversion. This is shown by the following: 

Assume there are n problem constraints. There are m on the first branch 
and m+1 through n on the second branch. We need examine only the diag- 
onal elements of A (because of symmetry). Let vector D 1 = diag (A) then, 
D corresponds element wise to the vector d'£' 



( 12 . 1 - 27 ) 


On the first branch of the trajectory the B matrix is defined by 



( 12 . 1 - 28 ) 
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and on the second branch, 

I" ^Savt 1 , m+1 


B = 


A 


n, n 


(12. 1-29) 


Clearly, on the main trunk of the trajectory (see Figure 11-1) 

B = A (12. 1-30) 

The matrix compression logic corresponding to the above description is pro- 
grammed in subroutine MTX3A. 

12. 2 PAYOFF IMPROVEMENT 

The concept of the expandable constraint vector is used to good advantage in 
PADS. During the first few iterations of a solution, only the constraint ele- 
ments are contained in the d^ vector. As soon as these elements are driven 
as close to zero as desired, the d'J' vector is expanded to contain the first 
improvement in payoff quantity, d$. Assuming this payoff or some fraction 
of it can be achieved on the next trial trajectory, a method is needed for 
estimating the payoff improvement for subsequent trials. 

As a result of asking for a d$>, a history of 6u and 6p changes result which 
will define the control metric (dP) 2 given in Equation (12. 1-11). Substituting 
Equations (12. 1-23) and (12. 1-24) into Equation (12. 1-11) and using the defi- 
nition of matrix A, 

(dP) 2 = (d¥.) T A _1 d¥. > (12.2-1) 
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It is now convenient to break A and ^ into minors. Let n be the number 
of constraints 


(d*.) = (d'k.'d$y r 


( 12 . 2 - 2 ) 



A l, n+1 


A 

n, n+1 
A , 

n+ 1 , n+ 1 


(12. 2-3) 


and 



(A 


n+ 1 , n 


n+ 1 , n+ 1 


We now rewrite Equation (12. 2-1). 


(12. 2-4) 


(dP) 2 


(d'k.'dtjj) A 



(12. 2-5) 


We can now expand Equation (12. 2-5) in terms of the partitioned matrix and 
vectors . 



A . . ,, d$ 2 + 2N T d^.’d$ + (d'k. T Md'k. ' 

n+1, n+1 i \ i / i 


( 12 . 2 - 6 ) 


This is a scalar quadratic equation in d<J> which has the solutions. 


d<$ 


-N T d*.' 

i 

A n+1, n+1 




n+1, n+1 


(12. 2-7) 
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It should be noted that the initial value of the metric, dP^, is calculated by 
the approximation 

,dP) 2 = d® 2 A n+1 , n+1 , < 12 ' 2 - 8 > 

This approximation assumes that the remaining elements of the d^ vector 
are all near zero. The equations above are programmed in subroutine 
PAY02. 

12. 3 ADJOINT EQUATIONS 

The adjoint differential equations are programmed in subroutine ADEQ3A. 
They appear in engineering notation in the same order as computed; that is, 
an ordered adjoint vector is one row in the adjoint variable matrix. Each 
row corresponds to one element in the (n+1) order d'I r vector. Hence, for 
the rth constraint (l<r<n+l), the adjoint variable or influence function for 
velocity, V, during a particular arc s is defined for convenience as 


x = x r S 

v v 


(12. 3-1) 


Subscript partial derivative notation will be used here. For example, 


Y = D * Y = 
v v dv 


(12. 3-2) 


It should be noted that the adjoint differential equations include only terms 
that can be non- zero. 


V X* + Y X* + h X* + m \ ‘ + + X* + P X* 

v v v Y vh vm v^ v p 


+ + 6 v X & ) 


X * 
Y 


= - (v X* + Y X* + h X* + 4* X* + P,, X* + fi v xu 
\ Y v Y Y Yh Y YP Y p/ 


(12. 3-3) 
(12. 3-4) 
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i* = -fy x* + y, x. 5 ' 5 + m, x.* + 4* x* + p, x'l + u, x* + 6, x*, 

h \ h v hY h m h h P ^h p hC 


(12. 3-5) 


=-(v x* + y x£ + m x.* + ^ x* ^ 

m \ m v "m^Y mm m ^ J 


(12. 3-6) 


4 


( V X.* + Y X.* + 4-, X* + P X* + u X*^ 

V 4> v 4 jY 4^+ ^ p 4< p/ 


(12. 3-7) 


(v„ X* + Y. x* + ^ Xf + u X- \ 
V p v py p 4< Pp/ 


(12. 3-8) 


X = 0 

P 


(12. 3-9) 


X Q . 0 


(12. 3-10) 


The partial derivatives of the equations of motion with respect to the states 
are not presented in detail here. They are programmed in subroutine 
PDY3A. The correspondence between engineering notation and FORTRAN 
symbol follows rules which facilitate interpretation. These are shown in 
Table 12. 3-1. 
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Table 12. 3-1 


EXAMPLES OF FORTRAN NOTATION 
FOR PARTIAL DERIVATIVES 


Engineering 

FORTRAN 

Engineering 

FORTRAN 

Symbol 

Symbol 

Symbol 

Symbol 

V 

VDV 

P 

0DV 

V 


V 



VDG 


UDR 

V, 

VDR 

6 

HTDV 

h 


V 


V 

VDM 


HTDR 


VDP 


VDA 

V P 

VD0 

K 

PDA 

V 

p 

VDU 


GDPH 

V 

V 

GDV 

•-3- 

PDPH 

h 

V 

HDV 

V 

1 a 

i V 

AVV > see Section 2^ 

m 

m 

MDM 

Y 

a h 

AGR see Section 2 j 

*P 

PDO 
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' Section 13 

STAGING OPTIMIZATION 

The rubber- stage design optimization equations which may be employed for 
Phase I sizing problems are described in this section. This type of optimiza- 
tion problem has a performance index (payload), which is a function of state 
values at more than one arc corner point, and as such is classed as a mixed 
boundary- condition problem. 


Equation (6.4-8) indicates that the payload is simply the difference between 
the final orbiter weight (burnout) and the empty structure weight. 

PL = W, - W (13-1) 

f e 

o o 

It is clear that is really the weight of the vehicle at final time 

o 



The orbiter structural weight is given in Equation (6. 2-2) as a function of 
orbiter propellant. 


W 


e 

o 



(13-3) 


We need to express the payload, PL, in terms of the mass (state vector 
element) at corner points where the variation of mass is non-zero and an 
explicit dependence of payload on mass exists. 


The weight of propellant used during the orbiter burn is 

w = w - w, 

P O f 

o 


(13-4) 
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At booster burnout, the booster empty weight is dropped, resulting in 


W 


W 


BO 


- W 


eB 


(13-5) 


where W „ is the booster burnout weight and W „ is the booster empty 
BO ^ ^ 

weight defined by 



W-n is the booster propellant weight, which in turn is related to the differ - 
”B 

ence between the vehicle gross lift-off weight and booster burnout weight. 



= W 


L.O. 


- W 


BO 


(13-7) 


Now, substitute Equation (13-7) into (13-6): 

W eB =f B (W L.O. ‘ W BO ) 

Then, substitute Equation (13-8) into (13-5): 

W o =W BO- f B ( W L.O. - W BO> 


(13-8) 


(13-9) 


(13-9) into (13-4): 


W p =W BO- f B ( W L.O. 

* r\ C 


(13-10) 


(13-10) into (13-3): 


W = f 
eo o 


W 


BO 


f B< W L.O. - W BO» 


W, 


(13-11) 
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and, finally, (13-11) into (13-1): 


PL = W £ 

o 



W 


BO 



(W 


L. O. 


W ) 
W BO' 



(13-12) 

(PAYLOD) 


This is what is needed to relate payload to the weight at the end of orbiter 
burn and at the end of boost. 

In order to construct the adjoint initial conditions or transver sality, condi- 
tions at the injection point, the explicit partials of payload with respect to 
mass are required. 


At orbiter injection 


8PL 
a m 


= (1 +f 0 ')g r 


(13-13) 

(PDBC) 

(PDBCQL) 


where 


f ' = 
o 



At booster burnout 


^L 
a m 


= - V <1 + f B '> g r 

BO 


(13-14) 

(ADJUMP) 


where 


9f B 

£ b’ aw 

p, 


Note that Q is invariant and f Q , f Q \ 

values of W * and W , respectively. 

P Q P B 


f B and f B * are 


evaluated at nominal 
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Section 14 

LAGRANGE MULTIPLIER TRANSFORMATION 


This analysis shows how the matrix of adjoint variables can be transformed 
into an approximation for the Euler -Lag range multipliers. These trans- 
formation equations are actually a useful by-product of proving that in the 
limit, the steepest descent solution satisfies the necessary conditions of the 
calculus of variations for the Mayer problem. 

We need to rephrase the problem of steepest descent in a manner similar to 
that presented by Denham and Bryson in Reference 8. We subdivide the con. 
straint vector and its corresponding sensitivities into d*. and d$and use as 

a goal the optimization of d<$ . We need to adjoin the expression for dP and 

the constraints. The resulting expression is [(see Equation (12. 1-12)]. 


T * r *04 

d$ = (s*) 6P + ^ / A J 6u j dt 

j-1 Jr i 


T 

+ v x d^ . 


) T W&Ujdt (14-1) 
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If this expression is regrouped 


d<$ = 


(S $ ) 


T, e *LT 
- v (S ) 


T 

M-(6p) Y 


6P 



T a 
v A 


.£ 2 . 

i J 


|j. (6u.) T wj 


6u.dt 

J 


+ v T d^ + H(dP) 2 


We take the first variation of d$> 


(14-2) 


6(d<J>) 


(S^) T - v T ( S ^ i ) T - 2fj.(6p) T Y 


6 (6p) 


i /( 

1 * / T . 


$ £2. „ <lr.£2. 

A J _ v A J - 2 hl(6 


j=l T 


Uj) T wj 6 (6Uj 


) dt 


+ v 6(d4'.) + (i. 6(dP) 2 


(14-3) 


2 

(dP) is a prespecified quantity; therefore, its variation will be zero. 
Likewise, d\lL is considered fixed and its variation is zero. For arbitrary 
variation of 6u. and 6p, Equation (14-3) can only equal zero if 



? ( S 


^i y 


- 2 (J.(6u.) Y = 0 
J 


(14-4) 


14-2 



and 


tffi. tf.fi. 

A J - vA J - 2h-(6u.) 


J 


W =■' 0 


(14-5) 


The parameter and control corrections are then 


(6p) 


= [<S*) T - V T 


(6u) = A J - w A 2 \i 


(14-6) 


(14-7) 


We may- 


rewrite Equation (14-7) as (see Equation 11- 1-16) 

\T 


T 

(6u.) = 


tffiT T 

-X J + v 1 \ X 3 ) 


D • f W" 1 
u 

2P 


(14-8) 


From 


the calculus of variations, we know the necessary condition that 


X D • f = 0 
u 


(14-9) 


Here, X is 


the Euler -Lagrange multiplier on the exact extremal. On the 


ion, 6uj is likewise identically zero. We may then infer that the linear 


solution 

transformation 


x?4 ,o i j-H*") 

is a good approximation of the Euler-Lagrange multiplier. 


(14-10) 
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By substituting Equations (14-6) and (14-7) into Equation (14-2) we may 
calculate In terms of the expressions used in Section 12. 2 [see Equa- 

tion ( 12. 2 - 3)] . 

v T = [n] T [m] (14-11) 


Henc e. 


x T .U* 0 j f- [N] T [M](^ in J) T (14 ' 12) 

where the adjoints are all functions of time and N and M are evaluated only at 
the beginning of the trajectory (at the end of the adjoint solution). This com- 
putation is performed in subroutine TRAN3. The stored adjoints are trans- 
formed and stored on sequential file for use later in the quasi-linearization 
portion of the program. 
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Section 15 

STEEPEST DESCENT NUMERICAL ANALYSIS TECHNIQUES 

This section of the report discusses some of the more important numerical 
analyses and convergence techniques employed in the steepest descent portion 
of PADS. Other aspects of the programming techniques used are found in 
Volume II. 

15. 1 NUMERICAL INTEGRATION, STORAGE RETRIEVAL, 

AND STEERING COMPUTATION 

The numerical integration method used in the steepest descent algorithm is 
standard fourth-order Runge Kutta, as described in Reference 9. The same 
method is used for both forward trajectory and adjoint differential equation 
integration. RKTA3A and RKTB3A are the integration subroutines. 

The basic integration cycle includes evaluation of the derivative twice at the 
mid-compute interval and twice at the end of the compute interval. During 
the integration of the adjoint differential equations, it is essential that data 
for the state vector be accurate at each evaluation of the adjoint derivatives. 
This is necessary because the adjoint derivatives are functions of the state 
on the nominal, f, trajectory. To accommodate this requirement, the 
adjoint integration has specially adjusted integration intervals to match up 
exactly with the time integration intervals of the nominal trajectory. An 
additional consideration arises at the mid-compute interval since the state 
is only an estimate there. To get a better representation of the state at the 
mid-compute interval, the following equation is used (ENTRY CORVAR in 
Subroutine REU3). 



At 


(15. 1-1) 


where y, is the state vector at the beginning of an integration interval, 
f |y Q | is the derivative at that point, and fjy Q + f(y Q ) ^-} is the derivative 
evaluated at an estimated mid-point state. 


During the forward trajectory integration, the following data are stored at 
each mid- and end-integration interval in preparation for the adjoint integra- 
tion: (1) y, the state vector; (2) a , <J> , the steering vector; (3) i, the arc 
number; (4) i p , the phase number; (5) K 1 flag, the option flag that tells which 
governing equation is to be used for thrust; (6) flag, the option flag that 
tells which governing equation is to be used for a; and (7) arc, phase, and 
elapsed time. 

During the integration of the adjoint differential equations, the forward tra- 
jectory data are retrieved at each mid- and full-interval to permit the 
calculation of the necessary partial derivatives. As the adjoint integration 
proceeds, a number of quantities are stored in preparation for computing 
the control and parameter corrections on the next trial trajectory. These 
include the adjoint matrix, \ 1 J; the impulse response matrix, A J ; and 

the A matrix excluding parameter contributions (see Equations (12. 1-21) 
and (12. 1-22)). 

At the mid-compute interval of the adjoint integration, no refinements of the 
adjoint variables are made. However, the impulse response matrix is 
refined before storing it. The reason for this selection will be given in the 
next few paragraphs. 

The optimized steering computation on the trial trajectory using retrieved 
data from the adjoint solution is a very critical aspect of the program. The 
technique that is used was developed in the course of trying a great number 
of different approaches and has proved to be the most stable and efficient. 

At the beginning of a compute interval at arc time t, the following informa- 
tion is available. 

A. y(t), the current state vector (trial trajectory) 

B. the old state vector (nominal trajectory) 
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• Q ■ 

C. \ 1 J(t), the adjoint matrix 

D. u P old' steering and parameter values (nominal 

trajectory 

E. S I, parameter sensitivities for parameters that have not yet been 
perturbed 

F. (t), the impulse response matrix. 

Equation (12. 1-23) is to be evaluated. The first calculation performed is 


6y (t) = y(t) - y old (t) 


(15. 1-2) 
(MTX) 


Then, if any parameters are still to be perturbed 
i / \T 

[ss]= s x y _1 vs j U 5 - 1_3 ) 

(MTX) 

At this point, the A and [SS] matrices are simultaneously added and com- 
pressed according to the extent of the d4q vector 


B(t) = A(t) - [SS] (15. 1-4) 

(MTX) 

Then B(t) is inverted using a symmetric matrix inversion method (Refer- 
ence 10). 


Using the compression logic again, the augmented constraint vector is 
calculated. 


RR(t) = dtf. 



(15. 1-5) 
(MTX) 


The first part of the calculation provides the TR vector. 


TR(t) = B(t) - 1 RR(t) 


(15. 1-6) 
(MTX) 
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At arc time t, the control may be calculated as 


u(t) = u (.t) + W -1 A* iQj (t) TR(t) (15. 1-7) 

(MTXt) 

It is necessary to calculate the control at the mid-interval and end-interval. 
The vector TR(t) is not recalculated. The equations are therefore 

A 1 ' A t~ — 1 At 

u (t + f) = U Qld (t + f*) + W A J(t + I 1 ) TR(t) 


and 


u(t + At) = u Qld (t + At) + + At) TR(t) (15. 1-8) 

As mentioned in Section 12. 1, the vector TR is recalculated only if the so- 
called closed-loop control mode is operating. Once the integration marches 
past the input elapsed time where the feedback stops, the recalculation of 
TR simultaneously stops and TR remains constant for the remainder of the 
trajectory. 

Optimized arc time corrections are calculated up until the arc time is actu- 
ally changed. At this point, the sensitivity is dropped from the matrix 
and its contribution to [SSj is eliminated. 


15.2 ARC CUT-OFF TECHNIQUE 

It is essential that the arc j cut-off function, £7 y he as numerically close to 
zero as possible and likewise that the state vector at tj be precise. 


When 

= t - Tj (15. 2-1) 

obviously the satisfaction of = 0 is trivial; however, when 

= G j( y (t))- (15.2-2) 
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where is the desired cut-off and G is an arbitrary nonlinear function of 
the relative states, some special considerations are necessary. 


At time t during arc j, one can estimate the value of G at the next integration 
interval 

Gj (t + At) = Gj (t) + G. (t) A t (15.2-3) 

If 

G.(t + At) > «*., (15. 2-4) 

J J 


then the cut-off condition is detected. It remains to get a good estimate of 
the integration interval required to satisfy 


G. (t + At') = n*. 
J J 


(15. 2-5) 


Let 


H = G (t) - G.(t - At) 

and 

H' = Q*. - G.(t) 

We may get a second order approximation of At'. Let 


* '(tf 

C = A2(A') 


(15. 2-6 


(15. 2-7) 


(15. 2-8) 


(15. 2-9) 
(15. 2-10) 
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B = A'(l. - A2) 

and 

A = 1. - B 
Then 

CH 

At' = A (t - At) + B(t) + t 

GjW 


(15. 2-11) 
(15. 2-12) 

(15. 2-13) 


This equation is programmed in subroutine DTF3. 

The integration will then march out an interval At 1 at which time it will be 
further refined. 


Let 

K = gJ t + At') - G(t) 

K ' = n*j - G j{‘ + At '} 

R, - K' 

R ■ K 

B' = (R 1 ) 2 (3 - 2R') 

C' = 1 - B' 

D' = R' (R* - l) 2 

and 

E' = (R') 2 (R> - 1) 


(15. 2-14) 

(15. 2-15) 

(15. 2-16) 

(15. 2-17) 
(15. 2-18) 
(15. 2-19) 

(15. 2-20) 
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Then, the refined time estimate for t is 


D 1 

t = C't + B'(t + At') + K 

j 


E' 

G. jt + At'} 


(15. 2-21) 


The entire state vector is refined through a similar third-order interpolation 
formula (Hermitian interpolation). This is programmed in subroutine 
YREF3. 

Experience with this algorithm has been very good with the exception of one 
type of situation. Suppose inequality [see Equation (15.2-4)] does not detect 
the cut-off. This would mean that at the next compute interval, the integra- 
tion would have to march backwards. This is not permitted owing to limita- 
tions of storage retrieval logic and therefore the cut-off equation will not be 
satisfied. The only remedy for this occurrence is to reduce the integration 
interval size to minimize the likelihood that it will happen. 

The fact that previous time-point data are used for this algorithm makes it 
imperative that more than two integration intervals be used in each arc. 

15. 3 SOLUTION CONVERGENCE LOGIC 

The basic sequence of convergence in the steepest descent program is 

A. Integrate first nominal trajectory 

B. Satisfy problem constraints 

C. Improve payoff only after problem constraints are satisfied 

D. Continue payoff improvement until predicted improvement is 
very small or iterations exceed maximum number. 

15. 3. 1 First Nominal Trajectory 

The first nominal trajectory uses an input control history (on cards) or a 
prior stored solution-control history (see Volume III). This trajectory is 
integrated and its state and control are stored as described in Section 15. 1. 
Upon completing the terminal arc cut-off, the nominal constraint misses are 
calculated and the constraint vector d4q is established. 


15-7 



15. 3. 2 Satisfaction of Constraints 

One or more iterations are normally required to satisfy the problem 
constraints. After the first nominal trajectory is completed (iteration 1) the 
first adjoint solution is integrated and stored. Using this stored data, the 
program will attempt to drive all of the constraint misses d*^ to zero. If 
these misses are relatively small to start with, they may be driven down to 
a smaller size or even less than tolerance levels. If, on the other hand, 
they are large, the control computation may diverge or the constraint misses 
may get even larger. In the event that divergence occurs, the asked-for 
constraint corrections will all be halved and another trajectory, using the 
same stored data, will be integrated. The prime criterion for a satisfactory 
constraint pass is 

|d * | < Tol (15. 3-1) 

(TEST) 


where Tol is a vector of input constraint tolerances or 


(15. 3. 2) 
(TEST) 


|d*. | £ Id*. | 

old 


This means that all elements of | d*^ J have decreased since the last nominal 
trajectory. If some have increased and others decreased, the final test is: 



(15. 3-3) 
(TEST) 


where n = number of constraints. 


This means that the sum square metric of relative constraint misses has 
decreased since the last trajectory. 

If after five halving trials these tests are not satisfied, the program will 
stop and go on to the next case. 

If, on the other hand, the test in Equation (15. 3-1) is satisfied, payoff 
improvement may begin. 
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15. 3. 3 Payoff Improvement 

Once the test in Equation (15. 3-1) is satisfied, the vector may be 
expanded to include d$ p the payoff improvement. 


The first payoff improvement is either an input quantity or can be estimated 
from an input of the expected final value of performance index <$. 

Subsequent program iterations will attempt to get further payoff improve- 
ment. If an asked-for d$ is very large, control divergence or payoff 
divergence may occur or constraints may be violated. In the event that any 
of these occur, d<£ will be scaled down by 1 /n/~ 2" and the metric (dP) will 
be halved. 

At each new payoff improvement iteration, the predicted d4> is calculated 
using Equation (12. 2-7). The new d<$ is compared with the old value to see 
if special convergence acceleration is appropriate. If so, the (dP) metric 
is scaled up to tend to give larger payoff improvements per iteration. This 
payoff improvement logic is programmed in subroutine PAY02. 

When the predicted d$ is smaller than an input minimum payoff improvement, 
the program will integrate the solution trajectory. (When a QE solution is 
flagged, the transformation of adjoints to Euler -Lagrange multipliers 
(Section 14) occurs before the solution trajectory integration. ) 
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Section 16 

CONTROL VECTOR SOLUTION - OPTIMAL CONTROL 

This section treats the determination of the optimal steering angles a and cj> 
in the QL trajectory module of PADS. 

Sections 16. 1 through 16. 5 establish the definitions of various quantities and 
notational conventions. In Section 16 . 6, the calculus of variations is applied 
to the trajectory optimization problem to derive the necessary conditions for 
optimality. Among these conditions are the optimal control laws for <J> and a 
which are discussed in greater detail in Sections 16. 7 and 16 . 8. Finally, in 
Section 16. 9 the effects of control and state variable inequality constraints 
are discussed. 

16. 1 THE STATE VECTOR IN THE QL MODULE OF PADS 

If the method of quasi-linearization were amenable to problems with variable 
end-points, the state vector in the QL module of PADS would be the same, 
except for order, as it is in the steepest descent module. However, to cir- 
cumvent the problemof variable end-points, QL makes use of a transforma- 
tion that incorporates an additional state variable. 

The transformation runs as follows: Let z denote the original state vector 

and suppose that z conforms to the differential constraint 

= g (t, z, U), t Q stst F (16. 1-1) 

where t is the independent variable (time), U is the control vector, tp is 
unknown, and g is a known vector function. Define 

t = t„ - t (16. 1-2) 

F o 
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and take 


x = (t - t Q )/T ( 16. 1-3) 

It follows that x(tp) = 1 and x(t Q ) = 0, regardless of the value of tp. Let us 
define a new state vector y* as 

y*(x) = z(t(x)), 0 < x < 1. (16.1-4) 


Then, by the chain rule of differentiation, 



dz dt 
dt dx 



0 < x < 1. 


(16. 1-5) 


Finally, we construct the actual state vector that QL treats by appending the 
variable t to the vector y*. Call this vector y. Then the quasitime deriva- 
tive of y is well defined as 


^ = f(x, y, U), 0 < x < 1 


( 16 . 1 - 6 ) 


where 


dy. dz. 

_i _ l 

dx dt 

if y. is a member of y* and 


(16. 1-7) 



(16. 1-8) 


i f Y, 


T 
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As a consequence, the state vector in QL is defined as 

y = (V, Y, iji, h, p, p, m, t, Q) T (16. 1-9) 

For a two-point problem — i. e. , a one arc problem — the variable t is simply 
the duration of the arc in seconds. In the next section, we will discuss the 
ramifications of multi-point problems. 

16.2 DESCRIPTION OF THE MULTI- POINT PROBLEM 

In the course of flying the trajectory, should anything happen to cause the 
right side of Equation (16. 1-1) to be discontinuous, then we have a multipoint 
problem. The point at which the discontinuity occurs is called a corner 
point. 

PADS recognizes four types of corner points: (1) branch points on the late 

side of which Equation ( 1 6 . 1-1) takes on more than one value; (2) end points 
which are the last points of branches; (3) initial points which are the first 
points of trajectories; and (4) intermediate points which are corner points 
that do not fall into any of the above categories. Since an intermediate point 
is not a branch point, the right side of Equation (16. 1-1) will have exactly 
one value on the late side of an intermediate point. 

If a problem has a branch point, it is called a branch problem. Presently, 
PADS can handle no more than one branch point in a trajectory, and the 
right side of Equation (16. 1-1) can have only two values on the late side of 
the branch point. 

The portions of the trajectory between the various corner points are termed 
subarcs. PADS assumes that all subarcs in a given trajectory are time- 
wise contiguous to some other subarc. In other words, there are no time 
gaps between any subarcs. 

Figure 16 . 2-1 depicts a typical branch problem in PADS. The subarcs are 
numbered sequentially so that for a branch problem the first subarc of the 
second branch is assigned the next number after the last subarc of the first 
branch. 
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As in the two-point problem discussed in Section 16. 1, the QL trajectory 
module transforms the time variable for the multipoint problem to the 
quasitime variable x. In this case, the transformation runs as follows: Let 

the corner points of the multipoint problem be ordered as follows: 

|‘n i + 1 < ‘n i+ 2 ^ (First Branch) 

V V *2 ••• *= ‘Nj < |t N2+1 <t N ^ +z ... < t Nj (Second Branch) 

( 16 . 2 - 1 ) 

For 0 + < x < 1 , define 

t (x) = t . - t (16.2-2) 

1 o 

for l + < x< 2 , define 

t (x) = t 2 - tj (16.2-3) 
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(16. 2-4) 


In general, for all I except I = N^-t-T and for (I- 1)"*" < x < I , define 


t (x) = tj - t ] ._ 1 


For I = N^ + l (branch problems only) and < x < (N 2 + I) , define 


T(x) = ‘n 2 +1 - ‘Nj 


(16. 2-5) 


Hence, for (1-1) < x < I , t(x) is the duration of the I tn subarc. 


16. 3 THE EQUATIONS OF MOTION 

In PADS, there are two different sets of equations of motion: the standard set 
which applies over most subarcs and a special set which applies only on 
vertical rise or pitch over subarcs. As a result, the equations of motion 
themselves are subarc-dependent, and, hence, for the multipoint problem, 
Equation (16. 1-6) becomes 

(16. 3-1) 
(NLDRV) 

where the subscript I on the right side of the equation indicates the subarc- 
dependency. 

16. 4 THE STATE BOUNDARY CONDITIONS 

A full description of the various boundary conditions in PADS that the state 
can be made to satisfy is given in Section 5. It is sufficient to note here that 
all of the state boundary conditions fail into two broad categories: state 

initial conditions and state target conditions. 


gj = fj- (x, y, U), (I-l) + < x< I" 


16-5 



Let if be the vector of all state boundary conditions for the entire multipoint 
problem. The QL module assumes that no state boundary condition involves 
the state at more than one corner point. Hence, we can partition if into 


if = 


;*o„; 


4^N 1 
'I'N +1 


^2__ 

4^n 2 +i 


c(;N 3 


| initial point 

intermediate points 

J branch point 

| intermediate points 

| end-point of first branch 

| intermediate points 
l end-point of second branch 


(16. 4-1) 
(MAGIC, 
BNDRY) 


The vector dj contains the state initial conditions at x = 0 + . 4 l contains the 

° + 1 
state initial conditions at x = 1 and the state target conditions at x = 1 . In 

general, at the intermediate point I, 4 * T contains the state initial conditions 

+ _ 

at I and the state target conditions at I . ^ contains the initial conditions 

_L t ^ __ 

at N 2 , the initial conditions at Nj and the target conditions at N^ - . 
contains the target conditions at N 2 " and contains the target conditions 

at No". 


Since the QL state vector y has nine components, then the vector if must have 
fewer than I 8 N 3 components; otherwise, no optimization could occur. 


16. 5 THE CONTROL VECTOR 

The control vector U in the QL module is defined as 

U = (T, 6 e> a, 4>) T (16.5-1) 
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There are several subvectors of U that have names of their own. The vector 


u = (a, (}>) T 


( 16 . 5-2) 


is termed the steering vector. On any given subarc, only the steering angles 
a and <j> may be optimized. The vector 

p = (T, 6 e ) T (16.5-3) 

is called the propulsion vector. Neither component of p may be optimized on 
any subarc. The vector 

w = (T, 6 e , a) T < 16 - 5 ‘ 4 > 


is called the in-plane control vector. As will be indicated later, the bank 
angle, , can always be determined independently of the vector w. This will 
be true both when is optimized and when it is not. As a consequence, the 
determination of the vector w will receive our greatest attention. 

In any case, in order to evaluate Equation (16. 3-1), it is necessary to deter- 
mine the vector U at every point of the trajectory. Let K be the vector of 
algebraic constraints which may totally or only partially determine the 
vector TJ, i. e. , 

K = K (x, y, U) = 0 06. 5-5) 

If K has four components, then U is completely determined and the vehicle is 
said to be undergoing non-optimal control. If K has less than four com- 
ponents, then U is only partially determined by K, and the vehicle is said to 
be undergoing optimal control. The minimum number of components K can 
have is two, in which case, both a and p are free to be optimized. 
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Clearly, both the number of components and the equations assigned to the 
components of K are subarc-dependent. Hence, Equation (16. 5-5) should be 
rewritten as 

K t = K t (x, y, U) = 0 (16. 5-6) 

1 1 (CONTRL) 

To completely determine U during optimal control, the indirect method of 
the calculus of variations is applied in the QL module of PADS. This will be 
discussed in detail in the following sections. 

16. 6 DERIVATION OF THE NECESSARY CONDITIONS FOR AN 
OPTIMAL SOLUTION OF THE MULTI- POINT PROBLEM 

Multi-point trajectory optimization to which PADS addresses itself is an 
example of the problem of Mayer. The problem statement runs as follows: 

We seek the state vector y(x) and control vector U(x) that cause the payoff 
function 

<t> = $ (y (No)) (16. 6-1) 

(BNDRY) 


to experience a minimum subject to the differential constraints 


dy 

as 



(x, 


y, U) 


0, 0 < x < N 3 


(16. 6-2) 

(NLDRV) 


and the algebraic constraints 


Kj. (x, y, U) = 0, 0 < x < N 3 


(16. 6-3) 
(CONTRL) 


and the fewer than 18N^ boundary conditions 


xjr = 


0 


( 16. 6-4) 
(MAGIC, 
BNDRY) 



Any valid target condition in the program can legally serve as a payoff. The 
fact that QL always minimizes the payoff should not concern the user, for if 
the maximization of the indicated payoff is desired, then QL minizes the 
additiv° inverse of the payoff. For example, suppose the maximization of 
m(N^) is desired. Then QL minimizes -m(N 3 ). 

There is one exceptional payoff function in the program that involves the state 
at more than the last point of the last branch. It is called the payload function 
and it is functionally dependent upon the booster and orbiter burnout masses; 
i. e. , 


$ = <t> (m (N “), m (N ")) (16.6-5) 

1 (BNDRY) 

This payoff applies only for Phase I sizing problems and will be treated 
separately. 

The approach to the problem stated above that the QL module of PADS 
employs is to adjoin the differential constraints and boundary conditions to the 
payoff to form the so-called augmented functional 


J = <E> + M • + 


n 3 r 1 

2 J x ' (ai - f i) ^ 


1=1 1-1 


( 16 . 6 - 6 ) 


where M and X. are vectors of Lagrange multipliers. In QL, the X. vector 
is defined as 


K = (\ 


V’ A -y» A -i r / 


m 


'■t' 


\q) 


(16. 6-7) 


and is referred to as the costate vector. 


Clearly, if y(x) and U(x) satisfy Equations (16. 6-2) through (16. 6-4), then 
Equation (16. 6-6) becomes 

J = $ (16.6-8) 
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Hence, the object of the QL module is to determine the y(x) and U(x) that 
satisfy Equations ( 1 6 . 6-2) through (16. 6-4) and cause Equation (16. 6-6) to 
have a minimum value. 


Before considering this problem, however, let us make use of the partition 
of SI/ introduced in Section 16. 4 in Equation (16. 6-6). We then get 


N 0 

N, r 

3 

3 r 

J = 4>+p • Ti + > p T * 

o ^o r I 

h + I j 

1=1 

1=1 1-1 


dx 


( 16 . 6 - 9 ) 


where the vector M has been partitioned into the subvectors p Q , p^, . . . , 
p^ 3 in the same manner as SI/. 

The first necessary condition for the minimization of J is that 

6J = 0 (16.6-10) 

where 6J represents the total variation of J. Expanding Equation (16. 6-10) 
we get 

6J = ^3<E>/ay(N 3 ')j 6y(N 3 ') + p Q • j^3 4 , Q /9y(0 + ) j 6y(0 + ) 


+ 1 iv [ 8 V 8 y< I+ > i 8 vw»] (l-yir-i) 

IeS ' \ ' 


+ ^N • 


9 4; N ^/9y(N 2 + )| 3 ^/SY^) | a+ N /ayfNp 


/6y(N 2 _ + j 
( 6y7N~+) 

XSy^-) 


+ H- n • 3^N 2 /3y(N 2 ') 6y(N 2 ") + p N • 3 4 | N 3 /3y(N 3 ' 
2r - 3 » 


6y(N •) 


n 3 .1 


+ 


6 / *•' (dx - £ l) dx 


( 16 . 6 - 11 ) 


L= 1 I- 1 
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where 


S = {I : 1 < I < N 3 and I / N , N 2 , N 3 } (16.6-12) 

The last term of Equation (16. 6-1 1) is the sum of the total variations of 

fch 

definite integrals. Let us consider only one of those variations — say, the I 
one. Furthermore, let us denote dy/dx by y'. Then, we have 


6 



(y' - fj) dx 



fj 6 y dx 

y 


(16. 6-13) 


Noting that 6(y') = (6y)', the first term on the right side of Equation (16. 6-13) 
can be integrated by parts to yield 


D 


I 


[ X * 6y] 


I 
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f T f y dx - / X. * f j 6 U dx 

y ii i u 

(16,6-14) 


where Dj denotes the expression on the left side of Equation (16. 6-13). 

Let the vector a consist of those steering angles that are free to be optimized 
on the ith subarc and let the vector q consist of the remainder of the U vector. 
Then 


U = (16.6-15) 

As pointed out in Section 16. 5, if Kj has four components, then o is the 
empty vector and q = U. If K^. has two components, then o = u and q = p. If 
Kj has three components, there are two possibilities depending on which 
steering angle is to be optimized. If <j> (the bank angle) is to.be optimized, 
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then cr = (cj>) and q = w. If a is to be optimized, then c r = (a) and 
q = (T, 6 e , <M T . 

In any case, because of Equation (16. 6-3), there is an implicit function k so 
that 


q = k (x, y, <r) 


( 16 . 6 - 1 6 ) 


(see Reference 11) and, as a consequence, Equation (16.6:- 14) can be 
rewritten as 


D 


I 


[\ • 6y] 


I 
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fj 6y dx 

y 




k 6y dx 

y 



fj k^ 6adx 


/' • ■ 
1-1 


6 c r dx 
(16. 6-17) 


After collecting the terms under the integrals in Equation (16. 6-17) on 6y 
and 6a, Dj becomes 



(16. 6-18) 


Since 6y and 6a appear under the integral sign, both of the integrals in Equa- 
tion (16. 6- 18) must vanish in order to satisfy Equation (16. 6- 10). Hence, by the 
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Fundamental Lemma of the calculus of variations (Reference 12) the Euler- 
Lagrange equations 



and 


- A. • 



0, (I- 1)+ < x < I- 


(16. 6-19) 
(NLDRV) 


(16. 6-20) 
(CONTRL) 


must be satisfied at 


each point of the I 


th 


subarc. 


Assuming that Equations (16. 6-19) and (16. 6-20) are satisfied., then Equa- 
tion (16. 6-18) reduces to 

D i = [ x ’ 6 y ]J _ 1 = . <16 \ 6 ‘ : 

and the last term on the right side of Equation (16. 6-11) becomes 



( 16 . 6 - 22 ) 


Still assuming that the Euler - Lagrange equations are satisfied on every sub- 
arc, Equation (16.6-10) in its expanded form becomes, after collecting terms 
at the corner points, 
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(16. 6-24) 
(COSTAO) 


where E . E,, E_, etc , are the row vectors 
old 

e o = h- 0 • [94> o /M 0+ )] - m° + ) T 


E 


I 


^ I /9y(I + ) | 94y9y(I-)j - X(I+) T | -X(I 


,IeS 


(16. 6-25) 
(COSTAI, 
INTRPT) 





9H* /9y(N, ) 

1 


a(n 2 +) 


X(N 1 + ) 


T ! 


I -X(N 


1 


(16. 6-26) 
(COSTAB, 
BRANPT) 




94* n /9y(N 2 -) 
2 


+ x(n 2 ") t 


(16. 6-27) 
(ENDPT) 




9 + N3 /9y(N 3 -) 


+ X(N 3 ") T + 


9 $/9y(N 3 -) 


(16. 6-28) 
(ENDPT) 


Equation (16. 6-23) is called the trans ver sality equation. Since the state 
variations at the corner points are independent, the trans verality equation 
will be satisfied only if E 0 , E^, . . . , Ej^ all vanish. Hence, the problem is 


to determine the values of jjl 0 , (j.j, . 
through (16. 6-28) to vanish. 


that cause Equations (16.6-24) 


Toward this end, we note that Equations (16. 6-24) through (16. 6-28) are of 
the general form 


p T A + b = 0 (16. 6-29) 
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where A is a matrix of partial derivatives and b is a row vector. For 
example, at the intermediate point I, p in Equation (l6. 6-29) is Pp 



and 


A 


94 . 1 /dy(I + ) 



(16. 6-31) 


The number of rows in p and A equals the number of boundary conditions that 
apply at the corner point. The number of columns in A and b is 9 for initial 
and end points, 18 for intermediate points, and 27 for branch points. Let us 
denote the number of rows and columns in A by m and n, respectively. Pre- 
sumably, the m rows of A are linearly independent row vectors. 

Should m = n, the solution to Equation (16. 6-29) is simply 


T 

H- 


-b A 


-1 


(16.6-32) 


However, when m < n, p is underdetermined, and, as a consequence, some 
additional boundary conditions must be added. These new conditions are 
called transver sality conditions and they involve both the state and the costate 


The ontogeny of the transver sality conditions runs as follows. Since the m 
rows of A are independent, they span a proper subspace of E n . Take B to be 
any set of n - m row vectors so that AuB spans the entire space E . This 
generates an invertible n by n matrix 


TA' 

A' - |“g- 


(16. 6-33) 


compute 



= -b 



(16. 6-34) 
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Then, if 6 = 0, it follows that 



T 

= p A 


T 

+ 0 B 



(16. 6-35) 


Hence, to satisfy Equation (16. 6-29) when men, we must determine y(x) and 
A(x) so that 9 , evaluated in Equation (16. 6-34) vanishes. 


In general, some trans ver sality conditions will exist at each corner point. 
So we define the vector © to be the trans ver sality conditions over the entire 
problem, i. e. 


0 


N 


1 


Nj + 1 


N 2 

0N 2 + 1 


N. 


(16. 6-36) 
(MAGIC, 
BNDRY) 


For those problems in which the special payoff (Equation (16. 6-5) applies, 
all we need do is add the partial of 4> with respect to m(N ^ ") to the right side 
of Equation (16. 6-26). 


To summarize the development to this point, recall that we began with a set 
of differential constraints [Equation (16. 6-2)] 


d£ 

dx 


fj(x, y, U) 


0 


(16. 6-37) 
(N LDRV) 
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to which we have had to add the differential constraints ^Equation (16. 6-19)] 



(16. 6-38) 
(NLDRV) 


We also had a set of algebraic constraints [Equation (16.6-3)] 


Kj (x, y, U) = 0 


( 16 . 6-39) 
(CONTRL) 


and to these we have added the algebraic constraints [Equation 


( 16 . 6 - 20 )] 


- A ■ 



0 


(16. 6-40) 
(CONTRL) 


Finally, to the boundary conditions [Equation (16.6-4)] 


'k = 0 


(16. 6-41) 
(MAGIC, 
BNDRY) 


we have had to add the trans ver sality conditions [Equation (16. 6-36)] 


0 = 0 


(16. 6-42) 
(MAGIC, 
BNDRY) 


The terms k y and in Equations (16. 6-38) and (16. 6-40) are, according to 
the implicit function theorem (Reference 11) 



and 


k 


<r 



(16. 6-43) 
(ALGCON) 


(16. 6-44) 
(ALGCON) 


16-17 



Thus, the calculus of variations has transformed our original optimization 
problem stated by Equations (16. 6-1) through (16. 6-4) into the multipoint 
boundary value problem stated by Equations (16. 6-37) through (16. 6-42). 

16. 7 THE OPTIMAL BANK ANGLE 

Substituting Equation (16. 6-44) into Equation (16. 6-40) yields 



If <t> (the bank angle) is a 
tion ( 16. 7- 1 ) is 




mponent of a, then the <t> 



(16. 7-1) 
(CONTRL) 


component of Equa- 


(16. 7-2) 
(CONTRL) 


However, as an examination of the various component candidates of in 
Section 10 will verify, when ^ is a component of a, none of the candidates is 
explicitly dependent upon ( i>. In other words, when <i>ea, = 0. Hence, 

Equation (16. 7-2) simplifies to 


-*• \ = - V'i - = ° (16.7-3) 

Plugging the expressions fory^ and 1 ^ 4 , into Equation (16. 7-3) and simplifying 
yields 

sin <J> - X'p cos <t> /cos y = 0 (16. 7-4) 


This equation has two solutions which are supplementary. Take 


C 




2 V 
cos Y 


(16. 7-5) 


Then 


sin <t> = \^/C 


(16. 7-6) 
(CONTRL) 
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cos <j> = A-y cosy/C 


(16.7-7) 

(CONTRL) 


and 

(16. 7-8) 
(CONTRL) 


( 16 . 7-9) 
(CONTRL) 

are solutions to Equation (16. 7-4). One solution will place <t> in the first or 
fourth quadrant and is therefore termed the belly-down solution. The other 
will place <t> in the second or third quadrant and is consequently called the 
belly-up solution. 

If both solutions are admissible, then we apply Pontryagin's maximum 
principle to guide us in the choice of the best solution (Reference 13). 
According to this principle, if the functional J in Equation (16. 6-9) is to 
experience a local minimum, then the Hamiltonian function as defined by 

(16. 7-10) 
(CONTRL) 

must experience a maximum. Hence, our choice of <2> will be the one which 
yields the largest value of H. 

One other important observation is to be made about 4>. The optimal bank 
angle is independent of the other components of U. As pointed out in Sec- 
tion 9. 1, the same is true of the non-optimal bank angles. Hence, on any 
subarc, we can solve for the bank angle explicitly, and the problem of deter- 
mining U therefore simplifies to determining the subvector w. 



cos <t> = -\y cosY/C 


sin <J> = -\f/C 
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16. 8 THE OPTIMAL ANGLE OF ATTACK 

If a (the angle of attack) is a component of cr, then the a-component of 
Equation (16. 7-1) is 



( 16 . 8 - 1 ) 

(A LI) 


In view of the observations made in the preceding section. Equation (16. 8-1) 
simplifies to 



( 16 . 8 - 2 ) 

(A LI) 


regardless of whether or not <t> is optimal. Of course, Kj should now be 
viewed as having only two components. 


At this point, we must consider what might happen to Equation (16. 8-2) when 
the total acceleration limit is encountered at some point in the I subarc 
while the vehicle is still in the atmosphere. Just prior to hitting the limit, 
Equation ( 1 6 . 8-2) is simply 


- 




K 


( 2 )' 

>E 


- 1 



= 0 


(16. 8-3) 
(AL1) 


where the subscript I has been dropped because it is understood and the 
superscript ^ indicates the second component of K. At the point where the 
limit is met, however, t 0, and Equation ( 1 6 . 8-2) becomes 



(16. 8-4) 
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Unfortunately, Equations (16. 8-3) and (16. 8-4) do not necessarily yield the 
same value of a at the point where the limit is met. A discontinuity of a, 
moreover, necessitates the introduction of another corner point. Worse yet, 
the value of T may exceed its maximum and, thereby require the imposition 
of the total acceleration limit as a control constraint on a. 


To avoid these difficulties, Equation (16. 8-3) is taken as the only form of the 
optimal control law for a. 

16. 9 CONTROL AND STATE VARIABLE INEQUALITY 
CONSTRAINTS 

There are three first-order state variable inequality constraints and three 
control constraints. In both cases, the angle of attack is chosen so that the 
constraint value is satisfied. Also, in both cases, the assumed underlying 
unconstrained control equation is the optimal control law [Equation (16. 8-3)] 


In the case of a state variable inequality constraint (SVIC), the constrained 
portion of flight must begin on the late side of a corner point because this 
portion of flight begins by matching the constraint as a boundary condition 
(Reference 14). Thereafter, the angle of attack will follow the SVIC until the 
optimal a resulting from Equation (16. 8-3) causes the time-rate of change . of 
the SVIC to be negative. For example, suppose the SVIC is the dynamic 
pressure, q. At the corner point, the constraint is matched; i. e. , q = q r 
Thereafter, if 


‘max 


^ (a ) > 0 
dt opt 


then a is chosen so that 


dq 

3F 


(a) = 0 


but if 


HF (a opt ) - 0 


then a is used, 
opt 


(16. 9-1) 


( 16 . 9 - 2 ) 


(16. 9-3) 
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Once off an SVIC, the program cannot go back on without the introduction 
of another corner point. 


Ordinary control constraints, on the other hand, do not require a corner 
point and, consequently, the constraint can go on and off at will. 

For both SVIC's and control constraints, at the point where the switch from 
constrained to optimal control occurs, it should be clear that Equa- 
tion (16. 6-37) is continuous. However, the continuity of Equation (16. 6-38) 
is not as obvious. In order to prove continuity, we first observe that at the 
switch point, both the optimal and the constrained angle of attack satisfy 


- X • 


f - f K" 1 
a P P 



0 


On the early side of the switch, Equation (16. 6-38) is 
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(16. 9-4) 


(16. 9-5) 


where C is the constraint. On the late side of the switch, Equation (16. 6-38) 
is 
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/f - f K -1 

1 y p p 



( 16 . 9 - 6 ) 


If Equations (16. 9-5) and (16. 9-6) are to be equal, then we must show that 
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Since K and C are nonsingular, the inversion in Equation (16. 9-7) is given 
by 
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Multiplying the inverse by the row vector X * 



f ) yields 
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Substituting Equation (16. 9-4) into Equation (16. 9-9) yields 
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Plugging Equation (16. 9-10) into Equation (16. 9-7) yields the desired result. 
Hence, both the state and costate come off the constraint surface tangentially. 
16. 10 SUMMARY 

In view of the results of Sections 9. 1 and 16. 7, regardless of whether the 
bank angle is optimal or non-optimal, we can always evaluate it explicitly. 

As a result, the problem of evaluating U reduces to solving for w. Section 9. 2 
discussed this problem when a is non-optimal. Section 16. 8 shows how 
Equation (16. 8-3) fills in for the missing component of Kj when a is optimal. 
Hence, we can view Kj as always having three components regardless of 
whether a is optimal or non-optimal. However, if we do, we also must view 
Kj as being explicitly dependent upon the co state; i. e. 

Kj. = Kj (x, y, A., w) = 


0 


( 16 . 10 - 1 ) 
(ALGCON) 



Section 17 


QUASILINEAR SOLUTION OF THE MULTI-POINT 
BOUNDARY VALUE PROBLEM 

Recall that the calculus of variations transformed the original optimization 
problem into a multi-point boundary value problem whose solution will 
indirectly yield a local minimum of the payoff function (Section 16. 6). The 
boundary value problem itself falls into the category of nonlinear first-order 
ordinary differential equations. 

The numerical method employed by the QL trajectory module of PADS is 
known as qua si -linearization. This is an iterative technique which is actu- 
ally an extension of the Newton-Raphson iteration to function spaces. Since 
a proof of the convergence of the iteration is beyond the scope of this docu- 
ment, the ensuing sections are intended as an exposition of how the method 
works in PADS rather than why. 

A simple two-point problem is discussed in Section 17.1. A number of 
observations are made about the nature of the quasi-linear solution and the 
multi-point problem is then addressed in Section 17. 2. The subsequent sec- 
tions discuss the significant mathematical and numerical problems and tech- 
niques that are attendant on the solution of the multi-point problem. 

17. 1 THE TWO-POINT PROBLEM 

Consider the system of first-order nonlinear ordinary differential equations 

(17. 1-1) 
(NLDRV) 

where 



dY 

dx 


= F(x, Y, w), 0 s x g 1 
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and where 


F T = (f (x,y, w) T - X- (f +fk)) (17.1-2) 

and where, according to Equation (16. 10-1), w is subject to the algebraic 
constraints 

K (x, Y, w) = 0 (17. 1-3) 

(ALGCON) 

together with the boundary conditions 

(* T i 0 T ) = 0 (17. 1-4) 

(BNDRY) 

Suppose Y is a solution to Equations (17. 1-1) through (17. 1-4), then provided 
F is twice continuously differentiable with respect to Y, Equation (17.1-1) can 
be written as a Taylor series 

= F (x, Z, W) + (Y - Z) + HOT (17.1-5) 

where Z = Z(x) is an element in the same function space as Y and W satisfies 
K (x, Z, W) = 0 (17. 1-6) 

and HOT denotes a second-order remainder term. 

The idea behind quasi-linearization is if Z can be chosen so that HOT is 
negligible, then the nonlinear system in Equation (17. 1-1) can be approxi- 
mated by the linear system 


9F 

9Y 


Y = Z 


ds 

dx 


F (x, Z, W) + 


(s - Z) 


(17. 1-7) 
(LINDRV) 



Because such a system is linear in s, it will have a solution of the form 


18 

s(x) = p(x) + ^ h i (x)c i 
i= 1 


(17. 1-8) 
(NOMINAL) 


where p(x) is some particular solution of Equation (17. 1-7) 



F (x, Z,W) + 



(p - Z) 


(17. 1-9) 
(LINDRV) 


hj (x), h^x), . . . , hjg(x) are a set of linearly independent solutions of the 
homogeneous differential equation 


dh. 
i 

dx 



(17. 1-10) 
(LINDRV) 


and Cj, c^, . . . , c^g is any set of scalars that causes the equation 


(*(s) T 1 0(s) T ) = 0 


(17. 1-11) 
(BNDRY) 


to be satisfied. 

Since HOT has been neglected, however, s will not generally be a solution 
to Equation (17. 1-1). As a result, an iterative process is employed wherein 
Z is replaced by s in Equation (17. 1-7) and a new s is computed. 

If this process converges; i. e. , if for any positive e, a positive integer M 
exists so that 


max 

0 £ x <■ 1 


s 

m 


(x) - s 


m+1 



< e 


(17. 1-12) 
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whenever m > M (m denotes iteration number), then it can be shown that 
it converges to a solution of Equation (17. 1-1). Moreover, it converges 
at a rate that is quadratic (Reference 15). 

Since the particular solution p(x) can be any solution of Equation (17. 1-7), 

we establish the following advantageous conventions in choosing the initial 

T • T 

value p(0). Some of the boundary conditions in ( ^ : 0 ) will be initial 

conditions on the state or costate. For example, 

m(0) - 50, 000 slugs = 0 (17.1-13) 

For each such condition, the state or costate variable involved is said to be 
fixed (known). All state or costate variables that are not fixed are said to 
be free (unknown). For all fixed variables, we set the appropriate compo- 
nents of p(0) to the known values. Thus, continuing with our example, the 
seventh component of p(0) is set to 50, 000. For all free variables, the 
appropriate components of p(0) will be set to the corresponding components 
of Z(0). For example, suppose A v (0) is free, then 

Pl 0 (0) = Z 1Q (0) (17. 1-14) 

Of course, on the first QL iteration, Z(0) is the value of the initial arc at 
x = 0 + . On subsequent iterations, however, as we have already noted 


Z(0) = s m-l (0) (17 - 1_15) 

The homogeneous solutions h^(x), l^x), . . ., h^g(x) must be independent. 
Hence, we use the following convention in choosing their initial values 


h.(0) 



1 if i=j 


0 if tfj 


(17. 1-16) 
(SALVE) 
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As a result of this convention, if H is the matrix whose columns are 
h.(x), i= 1, .... 18, then 


H (0) = I 


(17. 1-17) 


the identity matrix. Moreover, 

s(0) = p(0) + H(0)c = p(0) + c (17.1-18) 

Since the values of the scalars c y c 2> .... c lg are chosen so that Equa- 

tion (17. 1-11) holds and since, by convention, p(0) is chosen so that it 
satisfies the initial conditions in Equation (17. 1-4), it is clear from Equa- 
tion (17. 1-18) that c i is trivially zero if Y £ is fixed. Hence, we can throw 
out those homogeneous solutions that correspond to fixed states. 

Thus the following becomes the actual convention for the homogeneous solu- 
tions: For each free variable introduce a homogeneous solution tu(x) 

whose initial value is 

h. (0) = 6.. (17.1-19) 

i iJ 

In actuality, then, the matrix H(x) will have 18 rows and up to 18 columns, 
depending on the number of free variables, n. 

Equation (17. 1-18) also conveys the meaning of the remaining c's. They 
are the necessary perturbations to the initial values of the free variables to 
cause s(l) to satisfy the target conditions in Equation (17. 1-4). Differen- 
tia-ting Equation (17. 1-8) by a particular c^ indicates the meaning of the 
homogeneous solutions. 


9s(x) = , , v (17. 1-20) 

9c. i ( ’ 

i 


These solutions represent the sensitivity of s(x) to a unit perturbation at the 
initial point of the corresponding free variable. 


17-5 



In view of Equation (17. 1-15), it is clear that one necessary condition of 
convergence in the sense of Equation (17. 1-12) is that 

lim c = 0 (17. 1-21) 

0 

Hence, as an alternate definition of convergence, we use 
n 

J |c.| = 0 (17. 1-22) 

i=l 

17.2 THE MULTI-POINT PROBLEM 

For a multi-point problem Equations (17. 1-1) through (17. 1-3) become 


(17.2-1) 

(NLDRV) 

(17. 2-2) 


Kj (x, Y, w) = 0 
respectively. 

More important, the following conventions are adopted for the particular and 
homogeneous solutions. 

If a variable is continuous across a corner point; i. e. , if 

Y i (I+) ' Y i (I_) = 0 (17. 2.4) 

is one of the boundary conditions in Equation (17. 1-4), then the correspond- 
ing component of the particular solution also goes across the corner 
continuously, 


(17. 2-3) 
(ALGCON) 


~ = Fj(x, Y, w), I - l + < x < i; I = 1, 2, .... N, 


F i = ^ f i^ x ' y > w ) 


- \ • / f, + f T k > 
I I y 

y q 


and 
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(17.2-5) 

(SALVE) 


p. (I + ) = p. (I“) 


If the mass state variable experiences a discontinuity of known or computable 
magnitude 


m (1+) - m (I ) + 5000 slugs = 0 


(17. 2-6) 


for example, then the seventh component of the particular solution experi- 
ences the same discontinuity 


p 7 (I + ) = p 7 (l“) - 5000 slugs 


(17. 2-7) 


If a variable Y. is free on the late side of a corner point, then set 
J 

p.(I + ) = Z.(I + ) 


(17. 2-8) 
(SALVE) 


and introduce a new homogeneous solution ln(x) whose initial value at x - I 
is given by 


^( 1 +) 


6 . . 

ij 


(17. 2-9) 
(SALVE) 


This is perfectly legal as long as we remember that the scalar for this 
homogeneous solution cannot be perturbed to satisfy target conditions prior 
to x = I + . 


(Because of the last rule, the number of homogeneous solutions increases 
monotonically as we progress from one subarc to the next. ) 

If a variable is fixed on the late side of a corner point; i. e. if 

Y i < I+ > - (vXe”) = 0 (17.2-10) 

is a boundary condition in Equation (17. 1-4), then set 
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known 

value 


(17. 2-11) 
(SALVE) 


Pi (I 



and zero out the i^ 1 row of the matrix H(I + ). This is valid because the i^ 
variable at x = 1^ has become insensitive to perturbations of free variables 
prior to x = I + . 


On a branch problem, if 


Y i ( N 2 ) - Y i( N i) = 0 


(17. 2-12) 


is a component of Equation (17. 1-4), set 


Pi ( N 2 ) = Pi( N i) 


(17. 2-13) 
(SALVE) 


If 


m 


^ - m + 5000 slugs = 0 


(17. 2-14) 


is a component of Equation (17. 1-4), set 


p ? (n+ ) = p ? (N^ - 5000 slug: 


(17. 2-15) 


If the mass is distributed between the two branches; i. e. , if 


m 


(N+j+m (n{). m (n‘) = 0 


(17. 2-16) 


is a component of Equation (17. 1-4) set 


P 7 (n+) = m (n‘)- m ( N ;) 


(17. 2-17) 


If a costate is distributed between the branches; i. e. , if there is a trans- 
versality condition such as 


V 


( n 2) + x v( n 1)-*v( n ;) = 


o 


(17. 2-18) 
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o 


in Equation (17. 1-4), then distribute the corresponding component of p. For 
example, Equation (17. 2-18) would result in 


p io ( N i>) = p io ( N i) - p 10 ( N i) 


(17. 2-19) 
(SALVE) 


The distribution of a costate results from the state going across the branch 
point continuously to both branches. 


The above conventions are chosen because they force 


s(x) = p(x) + H(x)c (17.2-20) 

to automatically satisfy all of the initial conditions in Equation (17. 1-4). As 
a result, the determination of the c's is based solely on matching the target 
conditions in Equation (17. 1-4). 

17. 3 RECOGNIZING THE INITIAL CONDITIONS ON THE COSTATE 
Since the transver sality conditions 0 are derived numerically, the program 
must have some means of recognizing those transver sality conditions that 
represent initial conditions on the costate. 

Let us first consider the initial point x = 0 + . At this point. Equa- 
tion (16.6-34) becomes 



= A (0 + ) T A* 


1 


th. / T T \ 

Suppose is the j component of ( p. q : 9 ). Then 


6 


o. 

i 


A (0 


j 


(17. 3-1) 
(COSTAO) 


(17. 3-2) 
(COSTAO) 
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where 


- 1 


V 


.th 


-1 


is the i column of A .If 
J o 


„ - 1 


= 6 k< ■ < 


0 if k i ( 


1 if k = t 


(17. 3-3) 
(COSTAD) 


for some f, 1 £ f £ 9» then 9 is of the form 

o . 

l 


«o = h (0 + » 

1 


(17.3-4) 


and the transver sality condition is 
\ ( (0 + ) = 0 


(17. 3-5) 


In other words, A f is fixed at the initial point. 

In fact, because the state variables at the initial point can only be fixed or 
free, the same will be the case for the costate variables. 

Next, consider the intermediate point I. At this point Equation (16.6-34) 
becomes 



Mi + ) 


T . T" 

: - A (I") 



1 


Suppose 0 j is 
i 


the j 


th 


component 


of 



Then 


9 


I. 

l 



- A(I') 



(17. 3-6) 
( COSTAI) 


(17. 3-7) 
(COSTAI) 
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Let us denote the top nine entries of 

*j 

by S and the bottom nine by T. Then, if T = 0 and S = 6^ for some t , 

1 < t S9, 



6 = A f (I + ) (17.3-8) 

i 

represents the fixed initial condition on A.^. If S = 6^ for some i, 

1 < I i 9, and T = S, then 

e 1 = X J (I + ) - X|(I") (17.3-9) 

i 

represents the continuous initial condition on A^. 


Finally, consider the branch point N^. At this point. Equation (16. 6-34) 
becomes 



(17. 3-10) 
(COSTAB) 


Suppose 9 i is the j 


th 


component 




Then 


6 . 
l 




(17.3-11) 

(COSTAB) 


Let us denote the top nine entries of 


.- 1 


N, 
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by R, the middle nine by S, and the bottom nine by T. Then, if T = S = 0 


and R = 6 , „ for some f , Is St s 9, 
kf 




H) 


(17. 3-12) 


represents the fixed initial condition on 



If T = R = 0 and S = 6, , 

k f 


9 . = 

l 



(17. 3-13) 


represents the fixed initial condition on 
then 



If S = 0 and R = T = 6 , 

kf 


B. = 
i 



(17. 3-14) 


represents the continuity initial condition on between the trunk and the 

second branch. If R = 0 and S = T = 6, . then 

kf 


9 . = 
i 


\( N l)-h( N 'l) 


(17. 3-15) 


represents the continuity initial condition on \p between the trunk and the 


first branch. Finally, ifR = S=T=6^, then 


6. = \ 
i j 


(f4) + y( N l)- h( N 'i) 

represents the costate distribution initial condition on Yj . 


(17. 3-16) 


The fixed, continuous, and costate distribution conditions are the only 
transversality conditions recognized as initial conditions on the costate. All 
other transversality conditions are treated as costate target conditions. 


17.4 SOLVING FOR THE C'S 

Having integrated the particular and homogeneous solutions forward to the 
end of the trajectory, we must then determine the values of the scalars 
Cj, c^, ..., c^ that, cause the target conditions in Equation (17. 1- 11) to 
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be satisfied. The method employed is the well-known Newton- Raph son 
interation for systems of equations. 


Let us assume that Equation (17. 1-11) has been reduced to target con- 
ditions only. Further, let us take the value 0 as an initial guess for the 

vector c. Define 

Ac (i) - c (i+1) - c (i) (17.4-1) 

where the superscript indicates the iteration number, and 

c^) = initial guess = 0 (17.4-2) 


According to the method of Newton-Raphson, the increment in Equa- 
tion (17.4-1) is given by 



a*/a c I" 1 /*(c (l) )\ 

_' 9 "e/9e J c=c (i) \e~(c (i) ")/ 


(17.4-3) 

(NEWCS) 


Clearly, as the vector of target misses on the right side of Equation (17.4-3) 
approaches zero, so will the increment Ac^ provided, of course, that the 
matrix of partial derivatives is non- singular. 

To see how the target misses and partials in Equation (17.4-3) are actually 
evaluated, consider the target conditions ^ at the corner point x = I. Sup- 
pose that prior to the point x = I, a total of m^. homogeneous solutions has 
been introduced. Let 



(17.4-4) 


be the vector of multipliers for these homogeneous solutions. In general, 
of course, the vector 
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(17.4-5) 



of multipliers for all of the homogeneous solutions in the problem will have 
more components than c^. Indeed 


c C 
1 


C 2 C 


C I C. . . C c 


If we let 


(17.4-6) 


Hj(x) - 


h j(x) '• h 2 (x) ’• 


• h (x) 

• m j 


(17.4-7) 


then 


sd") = p(l") + Hjd") Cj (17.4-8) 

and, by the chain rule of differentiation, 

H I /ac I = (a^/asd")) osin/acj) 

= (8« t /8s(I“)) H (I“) (17.4-9) 

(INTRPT, BRANPT, ENDPT) 

Of course, the partial derivatives of 'fj- with respect to those c's in Equa- 
tion (17.4-5) that are not in Equation (17.4-4) are zero. 

Since the transver sality conditions at x = I are numerically derived, their 
partials with respect to the c's are numerically approximated by divided for- 
ward differences. For example, if c^ is one of the components of the vector 
Cj, then 

8<V 8c. ~ (9 (s(c +Ac ) - 0 (s(c )))/A c (17.4-10) 

(INTRPT, BRANPT, ENDPT) 

Of course, 

sfCjiAc.) = s( C j) + Ac.h.d') (17.4-11) 
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Some of the transver sality conditions in 9^ may involve both sides of the 
corner point. Hence, the partials with respect to the c 's introduced on the 
late side of the corner point must also be approximated by divided 
differences. 

17. 5 EVALUATION OF THE SYSTEM JACOBIAN 

In order to numerically integrate Equations (17. 1-9) and (17. 1-10), we 
must evaluate the system Jacobian 



Recall that 


V T . T i T. 
Y = (y i X ) 

i 


and 


( 


F T = ( f T ! - A • (f + f k ) 

i y q y 


) 


Let 


and 


X' 


T 


“y 


+ f k ) 

q y 


Then 


(17. 5-1) 


(17. 5-2) 


(17. 5-3) 


(17. 5-4) 


9F _ 9y'/ay [_9y7dX_ 

9Y jy=z Lax'/ ay i ax'7axj Y=z 


(17. 5-5) 
(NLDRV) 


■ Naturally, the first step in evaluating the Jacobian is to solve Equa- 
tion ( 17 . 1 -3 ) for w. This problem has already been discussed for nonopti- 
mal control modes in Section 10.2. Nevertheless, let us see what the four 
submatrices turn out to be when both a and <j> are non-optimal. 
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When Q is non- optimal, q = w. 


Hence 


ay ' = 

ay 


f + f k 

y w y 


+ fj 


<b 


y 


(17. 5-6) 
(NLDRY) 


As we noted in Section 10. 1, there are only two non-optimal bank angle 
modes. For the mode in which 

4> = 0 (17.5-7) 

the term <j) vanishes. For the vertical rise or pitchover mode, the term 
vanishes. Hence, Equation (17. 5-6) simplifies to 

-111 = f + f k 
ay y w y 

Since q = w, Equation (16.6-16) becomes 

w = k (x, y) (17. 5-9) 


(17. 5-8) 
(NLDRV) 


As a result, the term k in Equation (17. 5-8) can be viewed as 



According to Equation (16. 6-43) 


(17. 5-10) 


8w 

ay 



(K + K , (j> ) 
y 9 y 


(17. 5-11) 
(ALGCON) 


When a and 4> are non-optimal. Equation (17. 5-4) becomes 


. (f + f |2f) = - x • 

'v way' a v 


ay 


(17. 5-12) 
(NLDRV) 


^Recall that 4> is an explicit function of y. , 
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A number of shortcuts become apparent. First, compute 9y'/9y by means 
of Equations (17. 5-8) and (17. 5-11). Then compute A' by means of Equa- 
tion ( 17. 5- 12). Since w and d do not depend on A, 


Ji' = o 

<3 A 

and 


- - rii!' 

8 A Idy. 


Moreover, the matrix 8A'/8y is symmetric. Observe first that 


f ax' "! 

i 

8A i 

9 

/ 8y'\ 

£ 1 

. a = - 8 Y -A 


3y i 

dY j 

av, 

8yj9y i 


Now 

9y' _ f + f 8w 
“ayT y. way. 


and, consequently, 


8 y 

8Yj 8Y£ 



9 f Q r 2 

^i , w 8w _|_ £ 9 w 

9yj dy^ 8y. w9yj9y i 


= f 


y.y. 


+ f 


8 w 
w 9 v- 
7 J 


+ f. 


w. 


+ f 


w wcfy 


9w 

-■57. 


0 


8w 

3Y; 


+ f 


d w 

w 8y. 9y. 
J 1 


The right side of Equation (17. 5-17) is symmetric with respect to 
Hence 


(17. 5-13) 
(NLDRV) 


(17. 5-14) 
(NLDRV) 


(17. 5-15) 


(17. 5-16) 


(17. 5-17) 


y . and y. . 
1 J 



(17. 5-18) 


2 2 
8 y* = 9 y' 

3yj3y i 9 y j3y i 

which implies, by Equation (17. 5-15), that 


i 

\ — 


' ax'" 

l9y : 

ij 

. 9 y . 


(17. 5-19) 


So we need only evaluate the upper triangular portion of 3X'/3y. 


Of course, we still must evaluate the terms 9 w/3y^3y^ for j > i. This 
evaluation has already been discussed in Section 10. 3. 


If 4> is optimal and a is non-optimal. Equation (17. 5-6) does not reduce to 
Equation (17. 5-8). Of course, Equation (17. 5-4) still holds, but Equa- 
tions ( 17 . 5-12) through ( 17 . 5-14) and ( 17 . 5-19) no longer hold. In this 
case 


8 y' 
3A 




and 


(17. 5-20) 
(NLDRV) 


i 




+ f 


; w 


w 3y. 


6 .. 

D 







X 


(17. 5-21) 
(NLDRV) 


Of course, ax'/9y is no longer symmetric. 


ax i 
3y j 


f + f — + f , 

y i y j wy i dy j 4> y i 



+ 


f 9w -|- f 
f wy. + WW 9 y. 4>y 


J 


J 


9w , , 9 w 

9 V W dy dy 

1 J 1 


. X 


(17. 5-22) 
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If both <i> and a are optimal, then q = p. Let K* be the first two components 
of K in Equation (17. 1-3). In other words 


K* 



(17. 5-23) 


The third component of K is given by Equation (16.8-3). Equation ( 17. 5-4) 
becomes 




+ f p ) 

p *y 


where 



(17. 5-24) 
(NLDRV) 


(17. 5-25) 
(ALGCON) 


In this case, the matrix 9y'/9y is still given by 


9y 


f 

y 


+ f A- + 

W ay 


£ b ^y 


(17. 5-26) 
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as in Equation (17. 5-6), and the matrix dw/dy is still given by 


9w 

97 


- K'i(K y + 


<t> ) 
y 


(17. 5-27) 
(ALGCON) 


as in Equation ( 17 . 5- 1 1) . However, the matrix 9y' / 9 A is now given by 


9y ' = f * 

9X w 8A ® A 


where 


9 w 
9A 


- K" 1 (K 

w ' 


A + K <b ^ 


(17. 5-28) 
(NLDRV) 
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(17. 5-29) 
(ALGCON) 



Moreover, the matrix 3 A 1 /3A is given, element-by-element, by 





+ f 

P 



6 . . 




+ f 

pu 




• \ 


(17. 5-30) 
(NLDRV) 


where, of course, u = (a ,<£). The matrix 3V /By is rather complicated and, 
needless to say, asymmetric. It is given, term-by-term, by 


3 \ i 


ay; 


3f. 3f 
_L + _Pp 
3Yj BYj P Yi 



A 


where 



3 w 

y. w3y. 
J 


+ f , . + 

y^ ^Yj 



+ f 


3w 

pw 3y. 

h 




P 3Yj 


• A 


(17. 5-31) 
(NLDRV) 










+ K 


3 w 

pw 3y. 
j> 



K 


y.y. 

7 i 7 1 



3w 



(17. 5-32) 
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In case a is optimal and <t> is non-optimal (i. e. , <j> - 0), then Equa- 
tions (17. 5-24) through (17. 5-32) still apply, but the terms <f> y and <t> A vanish 


17. 6 METHODS OF INTEGRATION 

Another major step in the procedure described in Sections 17. 1 and 2 is the 
integration of Equations ( 17 . 1 - 9 ) and ( 17. 1- 10). This is accomplished by 
means of a fourth-order Adams-Moulton scheme with a standard fourth- 
order Runge-Kutta starting procedure. 


Let F(x, s) denote the right side of Equation (17. 1-9) or (17. 1-10). Then, 
for the Runge-Kutta starting procedure, we have 


s(x+h) 


s(x) + ^ + 2 k-, + 2 k^ + k^J 


(17. 6-1) 
(RKUTT1, RKUTT2) 


where h is the step size and 


kj = F[x, s(x)], 


(17. 6-2) 
(RKUTT1, RKUTT2) 


k 2 = F^x + |, s(x) + | k } j 

k 3 = F [ x + I’ s(x) + l k 2] 


(17. 6-3) 
(RKUTT1, RKUTT2) 


(17. 6-4) 
(RKUTT1, RKUTT2) 


k 


4 


F 


+ h, 


s(x) + h k^ 


] 


(17. 6-5) 
(RKUTT1, RKUTT2) 


This starting procedure is applied over the first three intervals of each 
subarc. The integration over the remainder of each subarc is accomplished 
by the Adams-Moulton process. This is a so-called predictor-corrector 
technique in which the predictor is given by 


s^x+h) = s(x) + 24 £^5kj " ^ k 2 + 37^3 - 9k^j 


(17. 6-6) 
(MADAMS) 
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where, once again, h is the step size, 
and 


the subscript 


P 


denotes predictor 


F[x, s(x)] 

(17. 6-7) 


(MADAMS) 

F C x - h, s(x - hfl 

(17. 6-8) 


(MADAMS) 

Fdx - 2h, s(x - 2h)J 

(17. 6-9) 


(MADAMS) 

F[x - 3h, s(x - 3h)] 

(17. 6-10) 


(MADAMS) 


The corrector is given by 


s c (x+h) = s(x) + p51i + 646k ^ - 264k 2 + 106k 3 - 19 k 4 j 

(17. 6-11) 
(MADAMS) 

where the subscript denotes corrector and 

c 

{ = F £x + h, Sp (x + h)J (Reference 9) (17. 6-12) 

The accuracy of the corrector can be improved by means of the following 
optional iteration. Let the superscript 1 denote the iteration number. Then, 
as the reader may verify, 

s c +1) (x+ h) = s^ l} (x+h) + zjz [251f (i) - 251f (i-1) ] (17.6-13) 

(MADAMS) 

where i= 1, 2, . . . , 


£ 


(i) 


x+h, s^ (x+h) 
c 


(17. 6-14) 
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and 


s<°> (x+h) = s (x+h) (17.6-15) 

c p 

We remark, in passing, that the error in the integration is proportional to 

A I „ 

h . Thus, if h s? 10" , the error will be on the order of 10 . On the other 

hand, the computation time varies inversely with the size of h. Hence, 
cutting the step size in half will double the run time. 

17.7 INTERPOLATION OF UNIVARIANT AND BI- VARIATE 
TABULAR FUNCTIONS 

In order to expand Equation (17. 1-1) into the Taylor series in Equa- 
tion (17. 1-5), Equation (17. 1-1) had to be twice continuously differentiable 
with respect to Y. Since Equation (17. 1-1) is functionally dependent upon 
both univariant and bivariate tabular functions, it is necessary to use 
interpolating functions for these tables that are twice continuously differenti- 
able. Two such interpolating functions are the cubic and bicubic spline 
functions. The former is used for univariant data and the latter for 
bivariate data. 

The cubic spline employed in PADS is the so-called natural cubic spline. 

Its derivation is given in Reference 16 and runs somewhat as follows. Let 
x be the independent variable and let u (not to be confused with the steering 
angles) be the dependent variable. Suppose the values of u are tabulated 
over the mesh 

A : x 0 < X 1 < ' • ’ < X N (17.7-1) 

Then, the table is given by the ordered pairs 

(x q , u q ), (xj, Uj), ..., (x-^j, u j^) (17.7-2) 

til 

Let M. denote the value of u"(x) at the i mesh point, i = 0, 1, . . . , N, and let 
1 

h. =x. - x. ,, i=l,...,N (17. 7-3) 

i i l-l 
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If we suppose that the second derivative is linear over each interval in the 
me sh, then 


u"(x) 


M 


i-1 


x. - x 
1 


h. 

i 


+ M. 
1 



X. . < X < X. 

1- 1 - - 1 


(17.7-4) 

(SPLINE) 


Integrating Equation (17.7-4) twice and evaluating the constants of integra- 
tion yields the results 


u(x) 


M. 


K - x ) 3 . ,, ( x - x iJ ; 


i-1 


6h. 


+ M. 


6h. 


+ 





(17.7-5) 

(SPLINE) 


and 


u'(x) 


•M. 


M 2 , w (*-* i-i) 


2h. 


+ M. 


2h. 


u. - u. . 
i i-i 

h. 


M. - M. , 
l i-i 


for x. , < x < 

l- 1 ~ 


x. . 

l 


(17.7-6) 

(SPLINE) 


The functions u(x), u'(x), and u"(x) defined by Equations (17.7-4) through 
(17.7-6) will be continuous at the mesh point x^ provided the quantities 
Mq, Mj, .... satisfy the relationships 


h. 

l 

6 


h. + h. , 

M. + — M. + 

i-i 3 l 


i+1 


M. 


i+1 


u 


i+1 

h. 


- u. 


i+1 


u. - 
1 


u 


i-1 


h. 


(17. 7-7) 


For i = 1, 2, ... , N- 1. At the end points we are free to close the values of 

M and M, T . We make the choice of M = M.. = 0. This implies a straight- 
o N o N 

line extrapolation outside the mesh. 
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Define 


h. 


i+1 


h. + h , 

1 l+l 


(17. 7-8) 
(MOMENT) 


n. = 1 - A. . 

r i l 


(17.7-9) 

(MOMENT) 


and 


d. 

i 


6 



h. + h 

i 


i+1 


(17. 7-10) 
(MOMENT) 


for i = 1, 2, . . . , N- 1. Then Equation (17.7-7) can be written as 

u.M. . + 2M. +\.M.. 1 = d. (17.7-11) 

r i i-I l i l+l l 

for i = 1, 2, . . . , N- 1 . 


Equation (17.7-11) together with the imposed end conditions on M q and M^ 
yield the linear system 


"1 

0 

0 

. . 0 

0 

0 


' M o ‘ 


‘0 

^1 

2 

1 • 

. . 0 

0 

0 


M 1 


d l 

0 

^2 

2 

. . 0 

0 

0 


M 2 


d 2 

0 

• 






• 


• 

0 

0 

0 

2 

X N-2 

0 


M N-2 


d N-2 

0 

0 

0 

^N-l 

2 

A N- 1 


m n-i 


d N-l 

0 

0 

0 

0 

0 

1 


L M n 


0 

m 





- 




- 


(17. 7-12) 
(MOMENT) 
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which is solved for M^, Mj, ■ . M.^. These quantities are then stored 

along with the table given in Equation (17. 7-2) so that the functions u(x), 
u'(x), and u"(x) given by Equations ( 17 . 7 - 5), (17.7-6), and (17.7-4), 
respectively, can be rapidly evaluated later on. 

The bicubic spline function employed in PADS is a natural extension of the 
cubic spline discussed above to functions of two independent variables. An 
excellent derivation of the bicubic spline interpolating polynomial has been 
given in Reference 17 and we will not attempt to duplicate it here. Suffice 
it to say that PADS' s and Reference 17 agree in every detail except for the 
following: 

Let the dependent variable u(x, y) be tabulated over the rectangular grid 


A x : x q < Xi < ... < x N 

A y ; y o < y l < ••• < y M 
In other words, we are given 

u„ = u (x^, y^ ), i = 0, 1, . . . , N; j = 0, 1, . . . , M 
Where Reference 17 uses 


(17.7-13) 


(17. 7-14) 


P = u x (x. y) 

q = u y (x, y) 


and 


s 


u 

xy 


(x, y) 


PADS uses 


p = “ xx <*.y> 


(17.7-15) 

(17.7-16) 

(17. 7-18) 
(17.7-19) 
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q = u (x, y) 

M yy 

(17.7-20) 

and 


s = u (x, y) 

xxyy 

(17.7-21) 

and while Reference 17 assumes one is given 


p. . = u (x., y.), i = 0, N; j = 0, 1, . . . , M 

X 1 'j 

(17.7-22) 

q.. = u (x., y.), i = 0, 1, . . . , N; j = 0, M 

y i i 

(17.7-23) 

and 


s.. = u (x., y.), i = 0, N; j = 0, M 

ij xy i 7 j 

(17.7-24) 

PADS arbitrarily takes ■ 


p.. = u (x., y.) = 0, i = 0, N; j = 0, 1, . . . , M 

F ij xx i' y y 

(17.7-25) 

(BLICO) 

q.. = u (x.,y.) = 0, i = 0 ,l,...,N;j = 0,M 

yy 1 J 

(17.7-26) 

(BLICO) 


and 

(17. 7-27) 
(BLICO) 

These are the bivariate analogs of the univariant end conditions - ( 

The quantities p.. at interior points of the grid are determined by univari- 
antly spline -fitting u.. along each grid line in the x direction. Similarly, 
the quantities q.. at interior points are determined by univariantly spline- 
fitting u.. along each grid line in the y direction. Finally, the quantities s^ 
ii J 


s.. = u (x.y.) 
ii xxyy Lj 


0, i = 0, N; j = 0, M 
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at interior points are determined by spline -fitting p.. along each grid line 
in the y direction. 


The coefficients of the bicubic polynomial in PADS are given by an equation 
that is similar to Equation 10 of Reference 17, but in PADS the A matrix is 
given by 


10 0 0 


A(h) 


- 1/h -h/3 1/h -h/6 

0 1/20 0 


0 - 1/ (6h) 0 l/(6h) 


(17.7-28) 

(BLICO) 
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